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ABSTRACT 

This  note  works  through  an  example  of  switching  between  many  coordinate  systems  using  a 
modern  matrix  language  that  lends  itself  to  describing  arenas  with  multiple  entities  such  as 
found  in  many  Defence  scenarios.  To  this  end,  it  describes  an  example  in  planetary  orbital 
theory,  whose  various  Sun-  and  Earth-centred  coordinate  systems  makes  that  theory  a  good 
test-bed  for  such  an  exposition  of  changing  coordinates.  In  particular,  we  predict  the  look 
direction  to  Jupiter  from  a  given  place  on  Earth  at  a  given  time,  highlighting  the  careful  book¬ 
keeping  that  is  required  along  the  way.  To  avoid  much  of  the  rather  antiquated  jargon  and 
notation  that  pervades  orbital  theory,  we  explain  the  first  principles  of  2-body  orbital  motion 
(Kepler’s  theory),  beginning  with  Newton’s  laws  and  proving  all  the  necessary  expressions. 
The  systematic  and  modern  approach  to  changing  coordinates  described  here  can  also  be 
applied  just  as  readily  in  contexts  such  as  a  Defence  aerospace  engagement,  which  follows  the 
interaction  of  multiple  entities  that  each  carry  their  own  coordinate  system. 
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Changing  Coordinates  in  the  Context  of  Orbital 

Mechanics 

Executive  Summary 

Real-world  defence  scenarios  might  be  described  or  managed  by  any  of  their  participants,  and 
a  core  element  of  this  description  is  the  ability  to  transform  between  the  many  coordinate 
systems  that  typically  quantify  the  entities  involved.  Switching  between  coordinates  is  often 
seen  as  a  classical  yet  difficult  problem.  This  report  attempts  to  show  that  the  task  can  be 
made  easier  and  more  transparent  by  using  unambiguous  notation  that  carefully  describes  all 
relationships  of  the  relevant  entities. 

A  worked  example  in  planetary  orbital  theory  is  a  useful  test-bed  for  such  an  exposition. 
The  various  Sun-  and  Earth-centred  coordinate  systems  involved  in  predicting,  say,  the  look 
direction  to  Jupiter  from  a  given  place  on  Earth  at  a  given  time  require  careful  book-keeping 
of  the  plethora  of  numbers  involved  in  the  calculation. 

With  that  worked  example  in  our  sights,  and  to  avoid  much  of  the  rather  antiquated  jargon 
and  notation  that  pervades  orbital  theory,  we  cover  the  first  principles  of  2-body  orbital  motion 
by  beginning  with  Newton’s  laws  and  proving  all  the  necessary  expressions.  The  main  focus 
here  is  to  show  how  to  interrelate  the  various  coordinate  systems  that  are  necessary  to  the 
worked  example. 

We  are  content  to  consider  the  2-body  problem  (Kepler’s  theory) — which  can  be  solved 
analytically,  unlike  the  many-body  problem — because  the  relevant  concepts  of  changing  co¬ 
ordinates  are  sufficiently  illustrated  in  a  2-body  scenario.  We  thus  decouple  the  Solar  System 
into  two  2-body  systems  that  are  gravitationally  independent:  Sun-Earth,  and  Sun-Jupiter. 
The  resulting  high  accuracy  in  the  prediction  of  Jupiter’s  look  direction  from  Earth  supports 
the  validity  of  this  decoupling. 

The  exposition  begins  with  the  relevant  classical  mechanics  and  time  concepts,  proves 
Kepler’s  three  laws,  then  establishes  and  describes  how  to  relate  the  different  coordinate 
systems  involved  with  the  Earth-centred  and  Sun-centred  inertial  frames,  the  Earth-centred 
Earth-fixed  frame,  and  the  observer’s  local  “flat  Earth”  frame.  It  describes  the  necessary 
celestial  geometry  and  orbital  elements,  and  finishes  with  the  worked  example  of  locating 
Jupiter  at  a  given  time. 

The  explanations  in  the  pages  that  follow  are  written  in  an  expansive  style  that  describes 
related  concepts,  such  as  what  equinoxes  and  solstices  are  and  when  they  occur,  how  julian 
days  are  defined,  and  how  orbital  calculations  can  be  extended  to  more  complex  motion,  such 
as  that  of  the  Moon.  In  particular,  the  theory  of  how  to  change  coordinates  more  generally 
in  an  elegant  but  also  powerful  way  is  explained  in  detail. 

Although  predicting  where  Jupiter  can  be  found  is  not  of  any  great  utility  in  a  Defence 
context,  anyone  who  understands  the  procedure  described  in  this  note  will  have  no  problem 
attempting  the  simpler  task  of  working,  for  example,  in  the  area  of  research  into  satellite 
positioning  systems  such  as  GPS. 
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1  Introduction 

This  note  serves  two  purposes.  First,  it  gives  a  worked  example  in  modern  mathematical 
language  of  the  subject  of  coordinate  transforms,  a  subject  that  is  the  acknowledged  bane  of 
many  a  researcher  in  Defence,  where  knowledge  is  sought  of  the  interaction  of  many  entities 
with  each  using  its  own  coordinate  system  to  communicate  its  view  of  a  scenario.  Second, 
this  report  describes  orbital  theory  in  a  streamlined  way  that  is  perhaps  more  in  line  with  a 
modern  physicist’s  way  of  viewing  the  subject  than  is  traditional  in  the  subject,  while  at  the 
same  time  avoiding  extraneous  concepts  such  as  the  theory  of  interplanetary  trajectories. 

The  world’s  increasing  reliance  on  satellite  technology  for  timing  and  positioning  suggests 
a  need  for  more  awareness  of  the  dynamics  of  Earth  satellites.  For  example,  in  the  study 
of  global  positioning  systems  such  knowledge  enables  a  researcher  to  construct  and  explore 
an  orbit  in  an  informed  way,  and  to  understand  the  relations  between  equivalent  choices  of 
parameters  that  describe  the  orbit. 

The  sections  that  follow  cover  the  relevant  classical  mechanics  and  time  concepts  from  the 
ground  up.  We  establish  the  necessary  frames  and  coordinates,  and  describe  how  to  relate 
the  different  coordinate  systems  that  are  required  to  calculate  a  satellite’s  position.  We  opt 
to  work  through  a  particularly  complicated  example:  finding  the  bearing  and  elevation  of  the 
planet  Jupiter  in  an  Adelaide  sky  at  a  given  time.  Although  Jupiter  is  hardly  a  satellite  of 
Earth,  the  method  of  predicting  its  position  covers  more  ground  than  that  of  forecasting  the 
position  of  an  Earth  satellite,  and  so  serves  as  a  more  in-depth  study  of  the  general  concepts. 
The  sections  that  follow  also  segue  into  descriptions  of  related  concepts,  such  as  details  of 
Earth’s  orbit  and  how  it  relates  to  the  seasons,  and  the  “julian  day”  approach  to  measuring 
time. 


2  Orbits  from  Newton’s  Laws 

Orbital  mechanics  is  an  old  and  established  subject.  Some  of  its  language  was  first  coined 
hundreds  of  years  ago,  and  the  subject  is  generally  still  presented  in  ways  that  reflect  its 
mediaeval  beginnings.  Its  literature  can  often  use  a  style  that  most  physicists  (myself  included) 
will  probably  regard  as  obscure:  for  example,  aside  from  sometimes  archaic  terminology, 
even  modern  books  might  treat  the  Sun  as  orbiting  Earth,  with  the  planets  orbiting  the 
Sun  [1],  This  choice  of  (non-inertial)  frame  is  useful  qualitatively  and  can  be  argued  for  and 
against  from  a  philosophical  point  of  view,  but  it  lacks  the  simplicity  of  an  inertial  frame  for 
calculating,  and  has  something  of  the  mediaeval  about  it.  It  is  rare  to  find  a  book  on  the 
subject  that  presents  succinctly  all  that  is  needed  to  locate  a  heavenly  body  without  being 
merely  a  “cook  book”  and  also  without  presenting  much  extraneous  information  at  the  same 
time  [2];  most  books  cover  a  lot  of  ground  over  many  chapters,  without  always  getting  to  the 
heart  of  the  subject  quickly.  The  treatment  of  vectors  and  coordinate  systems  in  many  books 
can  benefit  from  following  the  more  modern  approach  used  in  [3]  that  is  specifically  tailored 
to  cope  easily  with  multiple  coordinate  systems. 

The  question  we  will  answer  in  this  report  is  the  following:  if  you  gaze  up  at  the  sky  tonight 
at  8  p.m.,  where  must  you  look  to  find  a  given  planet?  Making  such  a  prediction  at  a  given 
time  and  place  involves  a  wealth  of  physics  and  mathematics,  and  seeing  the  calculations  come 
together  is  nothing  less  than  doing  an  experiment  in  classical  mechanics  on  a  grand  scale.  In 
this  report  I’ll  show  precisely  how  to  do  that  calculation,  by  building  the  necessary  theory 
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from  first  principles  and  then  joining  this  to  a  set  of  published  orbital  elements  describing  the 
planet’s  orbit  that  have  been  produced  from  astronomers’  observations. 

Predicting  where  a  planet  will  be  seen  at  some  time  and  place  follows  time-honoured 
classical  mechanics.  First  use  “force  equals  mass  times  acceleration”  to  analyse  a  situation 
governed  fully  by  gravity,  to  find  the  planet’s  position  in  its  orbit  as  a  function  of  time, 
using  appropriate  initial  conditions.  Orientate  this  orbit  (together  with  the  planet)  correctly 
relative  to  the  Sun,  do  the  same  for  Earth,  and  sight  the  planet  from  Earth.  Finally,  express 
the  resulting  vector  as  a  bearing  and  elevation  for  the  given  place  of  observation  on  Earth. 

Modern  fast  computers  calculate  planetary  ephemerides  (tables  of  predicted  positions  over 
time)  by  incrementally  solving  many  differential  equations  for  the  n-body  problem  that  com¬ 
bine  the  many  gravitational  influences  on  and  due  to  the  planets.  We  would  have  no  choice 
but  to  calculate  this  way  if  we  required  highly  accurate  predictions.  In  contrast,  treating  the 
chosen  planet-Sun  pair  as  a  2-body  problem  in  gravity  makes  it  exactly  solvable,  and  the 
result  is  sufficiently  simple  and  accurate  to  give  an  appreciation  for  the  physics  involved.  The 
resulting  keplerian  orbit  forms  the  basis  of  all  orbital  theory.  We  will  make  the  approxim¬ 
ation  that  for  the  purpose  of  predicting  Jupiter’s  direction  from  Earth,  the  Solar  System’s 
dynamics  can  be  decoupled  into  two  2-body  systems  that  are  gravitationally  independent: 
Sun-Earth,  and  Sun-Jupiter.  Each  of  these  2-body  problems  can  then  be  solved  individually. 
This  approximation  works  extremely  well  in  practice.  Additionally,  the  2-body  treatment  is 
completely  sufficient  to  serve  as  a  platform  for  discussing  the  necessary  coordinate  transforms. 

In  the  2-body  problem,  the  planet  is  subject  to  only  the  gravitational  force  exerted  by 
the  Sun.  It’s  an  undergraduate  task  in  mechanics  to  show  that  the  motions  of  the  centres  of 
mass  of  these  assumed-spherical  bodies  is  identical  to  the  motions  calculated  by  replacing  the 
planet  and  Sun  with  two  point  masses  that  are  located  at  the  centres  of  the  spheres  and  have 
those  spheres’  masses.  We’ll  assume  this  result. 

Now  consider  such  a  “point”  planet  of  mass  m  interacting  with  its  primary  (say,  the  Sun) 
of  mass  M ,  and  measure  the  planet’s  position  r  relative  to  its  primary.  The  “central  force” 
exerted  by  the  Sun  implies  that  the  planet  moves  in  a  plane,  and  so  we  place  the  Sun  at  the 
origin  of  orbit-plane  coordinates,  a  cartesian  set  x0P,  yOP,  zOP.  (These  are  more  conventionally 
called  perifocal,  or  PQW  coordinates).  Begin  by  analysing  the  planet’s  motion  in  the  orbit 
plane  xOPyOP  using  polar  coordinates  r,9,  where  r  =  |r|,  as  shown  in  Figure  1. 

Every  point  in  the  orbit  plane  has  attached  to  it  two  unit  vectors:  the  radial  unit-length 
vector  ur  points  everywhere  radially  outward  from  the  Sun,  and  the  transverse  unit-length 
vector  ud  is  produced  by  rotating  ur  by  90°  in  the  orbit  plane  right-handed  around  zOP. 
Referring  to  Figure  1,  we  apply  Newton’s  law  F  =  ma  to  the  mass-m  planet: 


—GMm 

- o - u 


position  of  m  relative 
to  any  point  fixed  in 
any  inertial  frame 


(2.1) 


where  G  is  Newton’s  gravitational  constant.  A  frame  is  defined  to  be  inertial  if  a  mass  released 
at  rest  in  that  frame  remains  motionless  indefinitely.  Any  frame  moving  at  constant  velocity 
relative  to  an  inertial  frame  will  itself  be  inertial.  A  convenient  choice  of  inertial  frame  for 
our  scenario  is  the  planet-primary  centre-of-mass  frame  (it  can  be  shown  that  centre-of-mass 
frames  are  inertial),  and  a  convenient  point  fixed  to  this  frame  is  the  centre  of  mass  itself. 
With  r  the  length  of  the  planet’s  position  vector  r  =  rur  extending  from  its  primary,  (2.1) 
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Figure  1:  Arrangement  of  vectors  and  masses  for  deriving  the  equation  of  motion  of  a  satellite 
of  mass  m  orbiting  a  primary  body  of  mass  M .  “CM”  denotes  the  centre  of  mass 
of  M  and  m:  it  divides  the  line  joining  the  masses  in  the  ratio  m  :  M .  Unit-length 
basis  vectors  at  three  representative  points  are  shown. 


becomes  (see  [4]  for  a  side  comment) 

—GMm  d2  Mr  Mm 

- - ur  =  m— 5-— - =  — - f  (2. 

r2  d t2  M  +  m  M  +  m 

(where  the  dots  denote  two  time  differentiations),  which  rearranges  to 

(2- 

where  /1  is  conventional  shorthand  for  G(M+m).  If  the  above  analysis  seems  straightforward, 
realise  that  we  have  derived  the  basic  equation  to  be  solved,  (2.3),  in  a  very  economical  way, 
due  to  our  combining  the  Sun-origin  coordinates  with  the  “CM-at-rest”  inertial  frame. 

We  will  solve  (2.3)  by  expressing  r  in  terms  of  ur  and  Ug,  and  then  equating  compon¬ 
ents  of  these  across  the  equals  sign  in  (2.3).  To  begin  doing  so,  we  establish  a  terminology 
that  will  become  important  later,  by  distinguishing  notationally  between  a  proper  vector  (an 
arrow  such  as  r),  and  its  description  in  some  coordinate  system  A,  being  an  array  of  co¬ 
ordinates  \t\a  known  as  a  coordinate  vector  [5].  So  write  r  in  the  orbit  plane’s  cartesian 
coordinates,  indicated  by  [r]0P: 


- G(M  +  m )  -u 

r  =  - ^ - ur  =  —FrUr 


r  cos  9 

rc 

r  sin  9 

rs 

(2.4) 


Differentiating  a  vector  is  easy  when  using  cartesian  coordinates:  just  differentiate  its  com¬ 
ponents  [6].  So  twice  differentiating  each  element  of  (2.4)  gives 


(r  —  rd2)c  —  (2  rO  +  rO)s 
(r  —  r62)s  +  (2  rO  +  rd)c 


(2.5) 
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Now  realise  that 

(2.6) 

and  combining  these  with  (2.5)  enables  r  to  be  expressed  without  using  cartesian  coordinates:1 

r  =  (r  —  r92)ur  +  (2 rQ  +  r6)ue  .  (2-7) 

Comparing  (2.3)  with  (2.7)  produces  the  two  key  equations  of  the  planet’s  motion  [8]: 

r  —  rd2  =  —  fi/r2  ,  2rQ  +  rd  =  0  .  (2-8) 

The  second  equation  in  (2.8)  can  be  rewritten  as 

=  (2-9) 

or  r20  =  a  constant,  conventionally  called  h.  Note  that  for  a  planet  of  mass  m  and  velocity  v 
relative  to  M,  the  angular  momentum  vector  per  unit  mass,  relative  to  M,  is 

L/m  =  rxv  =  rxr.  (2.10) 


Remember  that  r  =  rur,  and  just  as  we  found  r  above,  we  can  also  differentiate  each  element 
of  (2.4)  just  once,  and  again  use  (2.6)  to  produce 

r  =  fur  +  r6ug  .  (2-U) 


Substitute  these  expressions  for  r  and  r  into  (2.10),  and  apply  the  distributive  law  to  the 
cross  product  to  give 

L/m,  =  rur  x  ( fur  +  r0ue)  =  r20uZop  =  h,uZop ,  (2.12) 


where  w,op  is  the  unit  vector  normal  to  the  orbit  plane.  So  h  is  the  planet’s  orbital  angular 
momentum  per  unit  mass  with  respect  to  its  primary.  Remember  that  because  the  planet  is 
conventionally  chosen  to  orbit  in  the  direction  of  increasing  9,  the  angular  momentum  L  will 
be  parallel  to  uz^,  making  h  always  positive.  We’ll  use  the  fact  that  h  is  the  zOP  component 
ofrxr  ahead  in  (4.3). 

Now  solve  the  first  equation  in  (2.8)  by  a  change  of  variables  u  =  1/r,  eliminating  t  in  favour 
of  9.  To  do  this,  note  that  t  can  be  treated  as  a  function  of  9  alone,  in  which  case  the  chain  rule 
of  differentiation  combines  with  the  definition  of  h  to  give  d/dt  =  d9/dt  x  d/d 9  =  hu 2  d/d 9. 
So  we  calculate  r  by  applying  two  time  derivatives  of  hu 2  d/d 9  to  1/u.  The  result  is 


r  =  -h2u2  d2u/d92  . 


(2.13) 


The  first  equation  in  (2.8)  becomes 

d2u  ix 

d02  +u~  p  • 

This  linear  differential  equation  is  easily  solved  to  give  u 
stants  B  and  90,  and  from  this  we  can  write 


(2.14) 


B  cos (9  —  9q)  +  n/h2  for  con- 


_ 1 _ _ _ P _ 

B  cos(9  —  00)  +  ix/h2  1  +  ecos(0  —  90) 


(2.15) 


where  p  =  h? / p  is  called  the  orbit’s  semi-parameter  and  e  its  eccentricity.  The  angle  90  merely 
specifies  some  bulk  rotation  of  the  orbit,  so  we  will  set  it  to  zero  without  loss  of  generality. 

1See  [7]  for  an  alternative  way  of  producing  (2.7),  that  doesn’t  rely  on  cartesian  coordinates. 
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Figure  2:  Basic  angles  and  axes  describing  an  elliptical  orbit.  The  dotted  circles  aid  in  visu¬ 
alising  how  the  eccentric  anomaly  E  is  defined  on  page  7. 

3  Kepler’s  Three  Laws 

When  0  ^  e  <  1,  (2.15)  with  90  =  0  describes  an  ellipse  in  the  x0F  yOP  plane  symmetrical  about 
the  x0P  axis,  shown  in  Figure  2.  This  is  the  path  followed  by  a  planet  with  insufficient  energy 
to  escape  its  primary,  since  values  of  e  ^  1  give  parabolae  and  hyperbolae,  open  conic  sections 
that  describe  only  unbound  objects  such  as  comets. 

The  orbital  ellipse  has  semi-major  and  semi-minor  axes  lengths  a  and  b.  Its  centre  of 
symmetry  is  the  origin  of  the  XY  axes  in  Figure  2,  so  that  it  satisfies  X2/a 2  +  Y2/b 2  =  1.  The 
Sun  occupies  one  focus,  at  the  origin  of  the  x0P  yOP  axes  where  X  =  ae.  The  semi-parameter  p 
is  the  value  of  r  when  9  =  90°,  and  Pythagoras’s  theorem  produces  p  =  a(l  —  e2). 

Setting  90  to  zero  corresponds  to  measuring  9  from  the  planet’s  perifocus:  its  point  of 
closest  approach  to  the  focus.  (The  perifocus  is  also  called  perihelion  for  orbits  about  the 
Sun,  and  perigee  for  orbits  about  Earth.)  This  setting  of  90  to  zero  is  always  done,  in  which 
case  9  is  called  the  planet’s  true  anomaly.  Here,  the  word  “anomaly”  isn’t  meant  to  imply 
that  something  is  wrong;  perhaps  it  originally  meant  a  departure  from  perifocus.  One  theory 
suggests  that  the  language  of  orbital  mechanics  was  deliberately  obfuscated  by  early  navigators 
to  prevent  widespread  knowledge  of  how  to  navigate,  thus  discouraging  would-be  mutineers  [9]. 

We’ve  found  that  the  orbit  of  a  planet  is  an  ellipse  with  the  Sun  at  one  focus:  this  is 
Kepler’s  First  Law.  In  general  (2.15)  describes  a  conic  section — and  only  a  conic  section.  A 
sometimes- found  misconception  holds  that  if  momentarily  perturbed,  a  planet’s  orbit  would 
collapse  and  send  it  spiralling  in  to  the  Sun.  But  this  cannot  be;  a  spiral  path  is  not  a  conic 
section,  so  cannot  describe  planetary  motion.  No  matter  how  slowly  it  moves,  a  point  mass  m 
that  is  set  to  move  not  purely  radially  can  never  fall  in  to  collide  with  a  point  mass  M ; 
all  motion  of  bodies  under  gravity  is  a  part  of  an  unending  orbital  motion,  unless  the  non¬ 
pointlike  nature  of  the  bodies  gets  in  the  way.  When  you  throw  a  rock,  then  neglecting  air 
resistance,  from  the  moment  it  leaves  your  hand  until  it  hits  the  ground,  it  is  orbiting  Earth’s 
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centre  in  a  tight  ellipse.  The  end  of  the  ellipse  along  which  the  rock  flies  before  it  hits  the 
ground  is  modelled  very  accurately  as  a  parabola,  but  the  full  path  is  an  ellipse.  (Of  course, 
Earth  satellites  do  eventually  fall  to  the  ground  if  they  orbit  low  enough,  due  to  atmospheric 
drag  that  imposes  a  force  which  we  have  not  included  above.) 

Related  to  this  idea  of  an  orbiting  rock  is  the  term  “orbital  speed”;  for  example,  the  orbital 
speed  of  a  satellite  of  Earth  is  generally  given  as  8km/s.  How  does  this  relate  to  the  rock  a 
few  lines  up?  Although  a  thrown  rock  is  in  orbit,  its  elliptical  orbit  intersects  Earth’s  surface, 
so  it  will  quickly  hit  the  ground.  The  orbital  speed  of  a  rock  is  something  of  an  idealisation, 
being  the  (constant)  speed  that  the  rock  needs  in  order  to  orbit  Earth  in  the  smallest  circle 
that  never  touches  the  ground. 

Figure  2  includes  a  shaded  area  A  that  is  swept  out  by  the  Sun-to-planet  vector  r  beginning 
at  perihelion.  Simple  geometry  shows  that  the  infinitesimal  orbital  area  dA  swept  out  by  this 
vector  in  a  time  interval  dt  is  d A  =  r2d0 /2.  It  follows  that 

dA/dt  =  r26/2  =  h/2 ,  (3.1) 


which  is  a  constant.  So  the  planet’s  position  vector  sweeps  out  area  at  a  constant  rate:  this 
is  Kepler’s  Second  Law.  Now  integrate  (3.1)  over  one  period  T  of  the  planet,  noting  that  the 
area  then  swept  by  the  planet  is  the  area  nab  of  its  elliptical  orbit: 


nab 


(3.2) 


It  follows  that  T  =  2nab/h.  Squaring  this  equation  and  eliminating  b  and  h  by  using  b  =  aV  1  —  e2 
and  h  =  y/Jcp  =  \J —  e 2)  produces  Kepler’s  Third  Law: 


rjn  2 


47T2 

G(M  +  m) 


(3.3) 


Each  planet’s  mass  is  much  less  than  that  of  the  Sun,  so  (3.3)  becomes  T 2  oc  a3  with  a  constant 
of  proportionality  that  is  approximately  the  same  for  all  the  planets,  as  found  empirically  by 
Kepler.  We’ve  extracted  Kepler’s  three  laws  from  Newton’s  theory  of  gravity,  and  indeed 
this  was  the  major  early  success  of,  and  support  for,  Newton’s  work.  Kepler  is  often  seen 
as  having  lived  in  the  mediaeval  era  “BN”  (“Before  Newton”),  but  Kepler  himself  came  very 
close  to  deducing  the  law  of  gravitation  and  creating  the  subject  of  calculus,  and  it  was 
very  much  partly  on  his  shoulders  that  Newton  stood.  In  fact,  Newton  was  criticised  by  his 
contemporaries  for  not  giving  Kepler  his  due  in  this  regard  [10]. 


4  Introducing  a  Time  Dependence 


Equation  (3.1)  re-introduced  time  into  the  description  of  the  planet’s  orbit,  after  we  had 
temporarily  replaced  f  by  0  just  before  (2.13).  The  increase  in  swept  area  with  time  is 
traditionally  described  in  more  detail  by  introducing  a  new  angle.  Because  the  ellipse  satisfies 
X2/a 2  +  Y2/b 2  =  1,  this  new  angle  E  can  be  defined  such  that 


Y  r  sin  6  X  ae  +  r  cos  6 

sin  E  =  —  =  — - —  ,  cos  E  =  —  =  - 

6  6  a  a 


(4.1) 
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The  parameter  E  is  called  the  planet’s  eccentric  anomaly ,  and  is  the  angle  from  the  X  axis  to 
the  red  line  drawn  from  the  XY  origin  in  Figure  2.  A  few  lines  of  algebra  yield  r  =  a(l  —  e  cos  E), 
which  together  with  (4.1)  means  that  from  E  we  can  extract  the  planet’s  polar  coordin¬ 
ates  (r,9).  But  how  does  E  depend  on  time?  Recall  the  usual  relations  x0P  =  rcos0  and 
Hop  =  rsind]  then  with  E  expressed  in  radians,  it  follows  that 

x0P  =  a  cos  E  —  ae  ,  y0P  =  b  sin  E  , 

x0P  =  -aE  sin  E,  y0P  =  bE  cos  E.  (4.2) 

Now  remember  from  (2.12)  that  h  is  the  zOP  component  of  r  x  v\  also,  because  in  the  full 
three-dimensional  OP  coordinates  we  can  write 

[r  x  u]OP  =  (x0P,  yOP,  0)  x  (x0pj  Von  0) 

=  (0,  0,  xOPyOP  —  xOPyOP) ,  (4-3) 

it  follows  that  h  =  xOPyOP  —  x0p2/op-  By  substituting  from  (4.2)  into  this  last  expression  for  h. 
we  rewrite  h  in  terms  of  E\  then  equating  with  h  =  2nab/T  from  (3.2)  produces 

E  (1  —  ecos  E)  =  2n/T  ,  (4.4) 

which  can  be  written  as 

(1  —  e  cos  E)  d E  =  27r  d t/T . 

Integrating  this  equation  from  t  =  tperi  when  the  planet  is  at  its  perifocus  (where  9  and 
both  zero),  we  arrive  at  the  celebrated  Kepler’s  equation  for  the  time  evolution  of  E: 

E  —  esin  E  =  2n(t  —  tperi)/T  .  (4.6) 

(Remember,  E  must  be  expressed  in  radians  when  using  (4.6),  even  though  we  are  free  to 
write  it  as  a  number  of  degrees  when  not  using  that  equation.)  Kepler’s  equation  tells  us  how 
the  planet  moves  along  its  orbit.  The  quantity  E  —  esin E  clearly  increases  uniformly  with 
time  from  0  to  2n  over  one  period,  which  gives  rise  to  its  name  the  mean  anomaly  M : 

M  =  E  —  esin  A  .  (4.7) 

Kepler’s  equation  (4.6)  can  then  be  written  simply  as 

dM/dt  =  2i t/T  .  (4.8) 

All  three  anomalies — true  9,  eccentric  E,  and  mean  M — start  at  0  at  the  perifocus  and  reach  27r 
one  period  later;  but  they  advance  at  different  rates,  and  only  M  advances  uniformly,  which 
is  precisely  why  M  exists  and  why  it’s  useful.  Unlike  9  and  E ,  the  mean  anomaly  M  has  no 
geometric  interpretation;  being  simply  shorthand  for  E  —  esin  A,  it  is  not  “natively”  an  angle 
any  more  than  sin  A  is  an  angle.  But  because  M  increases  uniformly  from  0  to  27T  over  one 
period,  it  can  be  treated  as  the  angle  traced  by  a  point  orbiting  the  Sun  (or  indeed  any  other 
centre)  at  constant  speed  in  a  circle.  So  in  that  sense  the  mean  anomaly  fulfills  the  dream  of 
the  ancients:  to  reduce  a  planet’s  motion  to  constant  speed  in  one  circle!  [11] 

As  a  side  note  on  Kepler’s  Second  Law,  because  both  the  mean  anomaly  M  and  the  swept 
area  A  increase  at  constant  rates  as  the  planet  moves  from  its  perifocus  at  time  zero,  the 
relative  swept  area  A/ (nab)  must  equal  M/(2n). 


(4.5) 
E  are 
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As  an  example  of  solving  Kepler’s  equation,  we  calculate  how  far  advanced  a  planet  is 
from  perifocus  (i.e.  the  value  of  9,  the  true  anomaly)  after  a  quarter  period  for  an  orbit  of 
eccentricity  e  =  0.1.  The  answer  is  not  90° — that’s  the  mean  anomaly!  Instead,  begin  by 
solving  (4.6)  for  t  —  tperi  =  T / 4,  taking  care  to  use  radians  when  solving  Kepler’s  equation: 

E  —  0.1  sinA  =  7t/2  .  (4.9) 

For  ellipses  (e  <  1),  Kepler’s  equation  is  always  easily  solved  by  writing  it,  or  (4.9)  in  this 

CclS6,  clS 

E  =  0.1  sinA  +  7t/2  ,  (4-10) 

which  is  then  iterated  from  an  initial  guess  of  E  that  suffices  always  to  be  zero:  that  is, 
begin  by  setting  E  =  0  in  the  right-hand  side  of  (4.10)  and  repeatedly  recalculating  E  using 
the  same  equation,  inserting  the  latest  value  of  E  into  the  right-hand  side  on  each  iteration. 
Remember  to  work  in  radians!  After  just  4  iterations  here  the  value  of  E  converges  to  1.6703, 
or  about  95.7°.  Now  (4.1)  gives  sin0  and  cos 9: 

sin0  =  -\J  1  —  e2  sin E  ~  0.990-  ,  cos0  =  -(cos E  —  e)  ~  —0.199-  .  (4.11) 

J *  f  T  'f 

The  unknown  positive  number  a/r  cancels  if  we  calculate  tan0,  and  9  is  clearly  in  the  second 
quadrant.  It  follows  that  9  =  101.4°  [12].  To  summarise,  for  this  eccentricity  of  e  =  0.1,  after 
a  quarter  period  the  three  anomalies  are 

M  =  90°,  E  ~  95.7°,  9  ~  101.4°.  (4.12) 


As  an  aside,  we’ve  found  the  planet’s  position,  but  what  is  its  velocity?  This  velocity  v  =  r 
is  already  known  in  OP  coordinates  from  (4.2): 


Mop  =  0R)P,  2/op)  0)  =  A  (-a  sin  A,  b  cos  E,  0) , 

so  with  E  given  by  (4.4),  the  velocity’s  OP  coordinates  are 

r  1  2-7r(— asinE,  6  cos  A,  0) 

™  =  T(l-ecosA)  ' 


(4.13) 


(4.14) 


Another  useful  equation  results  if  we  form  the  dot  product  of  (2.3)  with  2 r.  An  integration 
and  some  manipulations  lead  to  the  vis-viva  equation  that  relates  the  planet’s  speed  v  (relative 
to  the  Sun)  to  its  current  distance  r  from  the  primary: 


v 


2 


(4.15) 


The  vis- viva  equation  applies  generally  to  elliptical  orbits  and  doesn’t  require  M  m.  In  fact, 
with  an  appropriate  definition  of  a  it  holds  for  all  conic-section  motion.  It’s  used  extensively 
in  the  art  of  manoeuvring  spacecraft  because  it  specifies  where  a  craft  should  switch  from  one 
orbit  to  another,  even  if  one  or  both  of  these  orbits  is  unbounded.  In  the  case  of  a  circular 
orbit  and  M  2>  m,  it  reduces  to  the  well-known  v2  =  GM/r. 
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4.1  The  Use  of  One  Global  Time  Variable 

Knowing  how  a  planet  moves,  we  wish  to  locate  it  at  some  given  time;  and  for  this,  a  single 
global  time  variable  t  is  a  necessity.  We  will  use  Greenwich  Mean  Time  (GMT)  [13],  but 
instead  of  using  a  24- hour  clock  that  resets  to  zero  each  midnight,  we  measure  t  in  days 
(including  any  fraction)  from  a  date  far  in  the  past  and  never  reset  it.  The  number  of  days 
elapsed  since  this  remote  starting  moment  is  known  as  the  julian  day  (JD) — no  relation  to 
the  julian  date ,  which  is  simply  a  date  in  the  julian  calendar.  The  JD  is  a  single  real-number 
way  of  expressing  GMT  and  so,  like  GMT,  defines  a  global  instant  [14]. 

The  remote  epoch  of  julian  day  zero  is  12:00  noon  GMT  on  an  arbitrary  date  chosen  far  in 
the  past  purely  so  as  to  render  julian  days  conveniently  positive  throughout  recorded  history. 
The  date  chosen  in  the  julian  scheme  was  1st  January  4713  BC  in  the  proleptic  julian  calendar , 
which  is  an  idealised  “well-behaved”  version  of  the  historical  julian  calendar — that  is,  without 
the  sporadic  changes  in  the  decreed  lengths  of  certain  months  that  occurred  in  the  historical 
julian  calendar  in  Roman  times.  This  start  moment  is  related  to  a  confluence  of  astronomical 
cycles  and  the  tax/census  cycle  of  ancient  Rome,  and  has  no  historical  significance.  In  the 
gregorian  calendar,  this  start  moment  was  12:00  noon  GMT  on  24th  November  4714  BC  [15]. 

The  JD  is  a  real  number,  so  typically  is  specified  to  several  decimal  places.  A  standard 
algorithm  to  convert  a  “year,  month,  day,  hour,  minute,  second”  to  a  julian  day  is  the  following, 
in  which  the  “floor”  function  returns  the  largest  integer  less  than  its  argument.  Set 

a  =  floor[(14  —  month)/ 12] , 
y  =  year  +  4800  —  a  , 

m  =  month  +  12a  —  3  .  (4.16) 

Then  with  UH,  M,  S”  =  hour,  minute,  second,  a  date  expressed  in  the  (proleptic)  julian  cal¬ 
endar  has  a  julian  day  of 

JD  =  day  +  floor  [(153  m  +  2)  /  5]  +  365y  +  floor(y/4) 

-  32,083  +  (H  -  12  +  M/60  +  S/3600)/24  .  (4.17) 

A  date  in  the  gregorian  calendar  has  a  julian  day  of 

JD  =  day  +  floor[(153m  +  2)  /  5]  +  365y  +  floor(?//4) 

-  floor(y/100)  +  floor(y/400)  -  32,045 

+  (H  -  12  + M/60  +  S/3600) /24.  (4.18) 

As  an  example,  apply  (4.16)  and  (4.18)  to  find  the  JD  for  the  gregorian  date  12:00  noon 
GMT,  1st  January  2000  (sometimes  referred  to  as  1.5  January);  you’ll  find  the  result  is  ex¬ 
actly  2,451,545.  You  can  verify  that  number  by  counting  the  days  from  first  principles,  being 
careful  to  remember  that  the  year  after  1  BC  is  AD  1:  so  4713  BC  is  AD  —4712,  a  leap  year. 

This  well-defined  time  variable  t  can  now  be  used  with  the  results  of  the  previous  sections 
to  compute  the  planet’s  motion.  An  epoch  t0  is  specified;  this  is  any  moment  that  astronomers 
agree  to  work  from,  usually  not  far  in  the  past.  The  position  of  the  planet  at  t0  has  been 
measured  by  astronomers  and  is  specified  as  part  of  the  set  of  parameters  that  describe  the 
planet’s  orbit.  We  ask  for  the  planet’s  position  at  time  f ;  the  answer  will  depend  wholly 
on  t  — 10.  An  easy  way  to  compute  the  time  t  —  t0  elapsed  since  the  epoch  is  to  convert  t0 
and  t  to  julian  days  JD(t0)  and  JD(f).  The  time  difference  t  —  t0  is  then  just  JD(f)  —  JD(f0) 
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Figure  3:  Relative  positioning  of  Sun- Centred  Inertial  coordinates  SCI  and  Earth-Centred  In¬ 
ertial  coordinates.  The  two  directions  marked  with  the  short  single  lines  are  the 
same,  as  are  the  pair  marked  with  the  short  double  lines.  Both  xsci  and  xECl  point 
to  the  First  Point  of  Aries  “T”  defined  in  Figure  4  and  in  the  text.  Earth  tilts  away 
from  zSCI  by  its  “polar  tilt” t  ~  23.439°.  The  small  blue  circle  next  to  zECI  shows  the 
direction  of  Earth ’s  precession  about  the  direction  of  zsci . 


“civilian  days”,  where  by  a  civilian  day  we  mean  normally  exactly  24  hours;  but  if  we  require 
one-second  accuracy,  all  leap  seconds  need  to  be  taken  into  account  too.  We’ll  ignore  leap 
seconds  in  this  report. 

Our  constantly  changing  galaxy  makes  it  impossible  to  predict  arbitrarily  far  into  the  future 
and  past,  of  course.  You  can  think  of  the  prediction  of  a  planet’s  position  at  time  t  as  a  zeroth- 
or  sometimes  first-order  Taylor  expansion  about  its  measured  position  at  t0.  We  will  assume 
that  most  parameters  don’t  change  with  time,  and  will  only  use  their  first  derivatives  when 
necessary. 


5  The  Earth- Centred  and  Sun- Centred 

Inertial  Frames 

Describing  planetary  motion  on  the  scale  of  the  Solar  System  requires  a  choice  of  frames  and 
coordinates  [16].  Two  useful  inertial  frames  quantified  by  cartesian  coordinates  are  shown  in 
Figure  3.  To  a  high  precision  their  coordinate  axes  are  unchanging,  pointing  to  almost  fixed 
points  on  the  celestial  sphere.  The  axes  of  the  Sun-Centred  Inertial  (SCI)  frame  xsci,  ysci,  zSCI 
originate  at  the  Sun.  Earth’s  orbit  defines  the  xsciysci  plane  (the  ecliptic  plane),  and  Earth 
orbits  right  handed  about  the  zSCI  axis. 

The  coordinate  axes  of  the  Earth- Centred  Inertial  (ECI)  frame  xECl,yECl,zEC1  originate  at 
Earth’s  centre.  The  zEGI  axis  is  Earth’s  spin  axis  pointing  out  of  our  North  Pole.  The  xEC1  axis 
is  parallel  to  the  xsci  axis,  and  the  yECl  axis  completes  the  right-handed  set.  Note  that  Figure  3 
deliberately  shows  Earth  with  no  detail,  to  emphasise  that  the  ECI  axes  bear  no  relation  to 
Earth’s  surface  features  such  as  its  countries. 

The  axes  of  the  Sun-Centred  Inertial  frame  are  shown  in  more  detail  in  Figure  4.  The  red- 
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Figure  4:  Earth’s  orbit  about  the  Sun,  showing  how  all  axes  are  defined.  Earth  carries  its  ECI 
axes  almost  without  changing  their  orientation  relative  to  the  stars,  apart  from  the 
slow  change  due  to  Earth ’s  precession  and  nutation. 


green-blue  axes  in  Figures  3  and  4  are  the  same  set:  ECI. 

Earth’s  polar  axis  zEC1  can  be  treated  as  a  vector  with  a  projection  onto  the  ecliptic  plane. 
This  projection  is  drawn  as  the  long  dashed  arrow  in  Figure  4;  its  length  has  been  increased 
for  clarity,  but  is  otherwise  irrelevant.  Throughout  the  year  this  projection  vector  makes 
an  angle  with  the  Sun-to-Earth  vector  that  varies  from  0°  to  360°.  The  angle  is  0°  in  late 
December  at  a  precise  moment  called  the  December  solstice,  when  southern  daylight  is  at  its 
longest  [17].  At  the  precise  moment  of  the  March  equinox,  again  by  definition,  the  angle  has 
increased  to  90°.  At  the  June  solstice  it  is  180°,  and  at  the  September  equinox  it  is  270°. 
These  key  moments  of  time  in  Earth’s  orbit  are  defined  by  Earth’s  tilt,  and  so  have  nothing 
to  do  with  the  orbit’s  major  and  minor  axes.  Earth  reaches  perihelion  in  early  January,  when 
the  Sun  is  5  million  km  closer  to  us  than  it  is  in  early  July.  The  Sun’s  closeness  in  January 
makes  its  apparent  area  7%  bigger  then,  than  in  July. 

These  four  instants  in  the  year,  two  solstices  and  two  equinoxes,  are  separated  by  90°  of 
true  anomaly  9;  they  are  not  separated  by  90°  of  mean  anomaly  M,  and  so  are  not  spaced 
exactly  3  months  apart.  But  the  eccentricity  of  Earth’s  orbit  is  low  enough  that  they  are 
spaced  almost  3  months  apart. 

The  xSCI  (equivalently,  xECI)  axis  is  defined  to  lie  along  the  Sun-to-Earth  vector  at  the 
moment  of  the  September  equinox.  This  axis  has  a  “vanishing  point”  in  the  sky  which  lies,  by 
construction,  on  the  celestial  equator,  being  the  projection  of  Earth’s  equatorial  plane  into  the 
sky.  This  vanishing  point  is  called  the  First  Point  of  Aries,  denoted  T,  the  most  important 
celestial  measurement  reference  used  by  astronomers  [18].  Why  Aries?  The  vanishing  point 
is  almost  fixed  in  the  sky,  but  precession  rotates  Earth’s  spin  axis  zECI  left  handed  about  zSCI 
once  every  25,770  years,  making  the  direction  of  a;SCI  (and  xECI)  change  by  360°  over  this 
period  [19].  The  First  Point  thus  moves  slowly  westward  through  the  zodiacal  constellations 
when  viewed  from  Earth.  In  around  2000  BC  the  First  Point  moved  from  Taurus  into  Aries, 
and  its  name  presumably  derives  from  astrology.  Precession  has  since  carried  it  into  Pisces, 
but  it  has  kept  its  old  name.  Over  the  next  few  centuries  the  First  Point  will  move  into 
Aquarius,  an  event  of  astrological  significance  that  spawned  the  catchy  1960s  pop  song  “The 
[Dawning  of  the]  Age  of  Aquarius”. 
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2/ecef 


Greenwich  sidereal  angle 


Figure  5:  Left:  axes  of  the  Earth- Centred  Earth-Fixed  frame  are  rigidly  fixed  to  Earth's  body 
in  precisely  the  positions  drawn.  The  xECEF  axis  is  located  at  the  longitude  of  Green¬ 
wich.  Right:  at  any  moment,  the  amount  by  which  Greenwich  has  turned  past  the 
First  Point  of  Aries  (xKCI)  is  the  “Greenwich  sidereal  angle”. 

6  The  ECEF  and  Sidereal  Angle 

The  xSC[  and  ysci  axes  slowly  rotate  as  Earth  precesses,  but  we  can  freeze  the  positions  of 
those  axes  at  some  epoch  and  retain  those  positions  for  all  future  times  of  interest.  When  we 
run  time  forward  or  backward  from  that  epoch  to  some  requested  moment  of  viewing  the  sky, 
these  epoch  SCI  axes  remain  fixed  by  definition.  Like  these  axes,  the  Earth-Centred  Inertial 
axes  x'KCI,  yKCl,  zKCI  remain  almost  fixed  in  space  over  time  spans  much  less  than  25,770  years. 
But  we  will  certainly  include  the  effect  of  precession  by  rotating  Earth  through  the  required 
angle  within  the  “frozen”  SCI  frame. 

Another  frame  of  use — this  one  non-inertial — is  the  Earth- Centred  Earth- Fixed  (ECEF) 
frame,  with  cartesian  coordinates  shown  in  Figure  5.  The  ECEF’s  axes  are  fixed  to  Earth;  the 
xECef  axis  points  from  Earth’s  centre  through  0°  latitude/0o  longitude,  the  zECEF  axis  coincides 
with  the  zBC1  axis  (Earth’s  spin  axis),  and  yECBF  completes  the  set.  Earth’s  spin  within  the 
Earth-Centred  Inertial  frame  manifests  as  the  ever-increasing  angle  between  the  “fixed”  xECI 
axis  (pointing  to  T)  and  “rotating”  xECEF  axis  (at  the  longitude  of  Greenwich).  This  angle  in 
Earth’s  equatorial  plane  is  the  Greenwich  sidereal  angle ,  and  increases  by  360°  in  the  time  it 
takes  Earth  to  rotate  once  with  respect  to  the  distant  stars,  being  one  sidereal  day ,  a  period 
of  23  hours,  56  minutes,  4.09890  seconds  [20]. 

Similar  to  the  Greenwich  sidereal  angle,  at  any  given  moment  the  local  sidereal  angle  of,  say, 
Adelaide  is  the  angle  in  Earth’s  equatorial  plane  from  the  xE0I  axis  to  Adelaide’s  meridian  (its 
great  circle  of  longitude):  this  is  just  the  Greenwich  sidereal  angle  plus  Adelaide’s  longitude. 
Knowing  the  current  local  sidereal  angle  equates  to  knowing  our  current  orientation  relative 
to  the  fixed  stars,  which  allows  us  to  determine  where  planets  and  stars  will  be  seen  in  our 
sky  right  now. 

To  find  the  current  local  sidereal  angle  we  need  a  start  point;  this  is  always  the  Green¬ 
wich  sidereal  angle  at  some  epoch  for  which  that  angle  has  been  measured  to  high  preci- 
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sion.  A  standard  epoch  used  here  is  “J2000.0”  which  uses  the  date  mentioned  in  Section  4.1: 
12:00  GMT  1st  January  2000,  at  which  moment  the  Greenwich  sidereal  angle — the  angle  in 
Earth’s  equatorial  plane  through  which  Earth’s  zero  meridian  had  turned  past  the  direction 
to  the  First  Point  of  Aries — happened  to  be  280.46062°.  (Remember  that  the  Greenwich 
sidereal  angle  at  time  t  is  a  unique  number  independent  of  epoch;  but  to  calculate  it  we  must 
start  somewhere:  we  must  surely  begin  with  an  epoch  for  which  we  know  the  angle’s  value 
accurately.) 

As  an  example,  what  is  the  sidereal  angle  of  Adelaide  at  local  (daylight  savings)  time 
9:00  p.m.  on  22nd  March  2014?  We  simply  add  three  angles  modulo  360°.  Begin  at  the  J2000.0 
epoch  when  the  Greenwich  angle  was  280.46°,  rotate  Earth  (and  hence  Greenwich)  from 
the  epoch  to  the  current  time  through  the  second  angle,  and  then  add  Adelaide’s  longitude 
of  138.60°.  The  time  through  which  we  must  rotate  Earth  is  the  requested  time  minus  the 
epoch.  Work  in  GMT:  convert  Adelaide’s  time  to  10:30  a.m.  GMT  22nd  March  2014,  then 
use  (4.16)  and  (4.18)  to  find  that  this  is  JD  2,456,738.9375.  The  time  elapsed  since  the 
J2000.0  epoch  is  this  JD  minus  2,451,545,  or  5193.9375  days.  The  angle  turned  through  by 
Earth  relative  to  the  First  Point  (modulo  360°)  since  the  epoch  is  then  its  hourly  rate  times 
the  number  of  hours  elapsed: 


360° 

23  h  56  m  4.09890  s 


x  5193.9375  x  24 h  ~  56.71°. 


So  the  sidereal  angle  of  Adelaide  at  local  time  9:00  p.m.  on  22nd  March  2014  is 


280.46°  +  56.71°  +  138.60°  =  115.77°  (modulo  360°  of  course). 

angle  at  epoch  Earth  turned  Adelaide’s  longitude 
_ i 

Greenwich  sidereal  angle 


(6.1) 


(6.2) 


At  this  time  the  First  Point  of  Aries  lies  115.77°  west  of  the  Adelaide  meridian,  measured  along 
the  celestial  equator.  Similarly,  at  any  time  and  place  on  Earth  the  local  sidereal  angle  is  really 
just  a  measure  of  where  we  see  the  stars  to  be.  As  the  stars  wheel  about  the  celestial  poles, 
any  one  of  them  can  be  visualised  as  being  at  the  end  of  a  giant  clock  hand  that  turns  through 
360°  in  the  above  23-56-4.09890  hours.  This  is  probably  what  has  prompted  astronomers  to 
call  the  sidereal  angle  the  sidereal  time.  With  this  name,  the  angle  is  usually  specified  in 
hours/minutes/seconds,  where  24  such  angular  hours  are  defined  to  be  exactly  360°;  that 
is,  24  sidereal  hours  of  angle  are  turned  through  by  Earth  relative  to  the  inertial  stars  in 
the  above  23-56-4.09890  civil  time  hours.  My  view  is  that  the  name  “sidereal  time”  is  an 
unfortunate  choice  for  something  that  is  simply  an  angle  and  is  always  used  geometrically; 
the  fact  that  the  spinning  Earth  acts  like  a  giant  clock  does  not  call  for  a  new  unit  of  time. 
Also,  the  use  of  a  time  unit  to  measure  this  angle  is  unfortunate  because  a  sidereal  hour  has 
no  useful  relation  to  a  civil  hour  of  3600  seconds,  and  even  when  “sidereal  time”  is  quoted 
in  hours/minutes/seconds,  it  still  denotes  an  angle,  not  a  real  time.  In  this  report,  “time” 
denotes  real  time  and  “angle”  denotes  angle,  and  all  hours  are  the  only  sort  that  I  maintain 
should  exist:  civil  time  hours,  the  sort  that  clocks  and  wrist  watches  measure  [21], 

On  a  related  note,  if  you  prefer  using  degrees  instead  of  radians,  do  write  your  angles  as  one 
decimal  number  of  degrees.  Just  like  pounds-shillings-pence,  expressing  angles  in  base  60  as 
degrees-minutes-seconds  is  an  old  but  clunky  practice,  one  that  only  obfuscates  calculations. 
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Figure  6:  The  elements  describing  an  orbit’s  orientation  relative  to  a  reference  plane  such  as 
the  ecliptic 


7  Orbital  Elements 

Having  located  a  planet  within  its  orbit  at  some  requested  time,  we  must  now  place  the  orbit 
in  its  correct  orientation  within  the  Solar  System.  The  orbit  is  typically  described  by  three 
angles  shown  in  Figure  6  that  refer  it  to  some  convenient  plane  such  as  the  ecliptic.  The 
point  where  the  planet  crosses  the  ecliptic  plane  in  the  direction  of  increasing  zSCI  is  called 
its  ascending  node,  denoted  Q>.  Two  angles,  the  orbital  inclination  i  and  the  longitude  of 
the  ascending  node  0  locate  the  orbit  plane  relative  to  the  ecliptic,  and  a  third  angle,  the 
argument  of  the  perifocus  oj,  locates  the  perifocus.  A  fourth  angle  (one  of  the  “anomalies”, 
usually  M)  specifies  the  planet’s  position  in  its  orbit  at  some  epoch  t0,  which  need  not  be  the 
same  as  the  epoch  used  above  for  calculating  sidereal  angle.  Together  with  the  semi-major 
axis  length  a  and  eccentricity  e,  these  angles  i,fl,oj,  and  M0  =  M(t0 )  comprise  a  set  of  six 
parameters  that  fully  describe  an  orbit  and  where  the  planet  is  located  in  it  at  the  epoch  t0. 
They  are  determined  from  observations  by  astronomers,  and  we’ll  take  them  as  given.  Their 
first  derivatives  are  also  available,  so  that  whereas  we’ll  obtain  enough  precision  by  taking 
e.g.  D  to  be  constant  for  each  planet,  we  could  also  use  fi(t)  —  fl(t0)  +  Cl(t0)(t  —  t0 )  for  higher 
precision.  Remember  that  the  longitude  H  for  planets  orbiting  the  Sun  is  an  angle  in  the 
ecliptic  plane  xsciysci,  measured  east  of  the  First  Point  of  Aries — as  opposed  to  the  everyday 
terrestrial  longitude  of  world  cities,  which  is  measured  east  of  the  Greenwich  meridian  in 
Earth’s  equatorial  plane  ^ecef^/ecef- 

Note  the  three  words  that  all  denote  an  angle  here:  longitude,  argument,  anomaly.  Perhaps 
these  words  date  back  to  those  maritime  practices  of  old  referred  to  earlier.  And  perhaps 
it’s  those  same  die-hard  practices  that  have  caused  oj  and  M0  to  be  sometimes  scrambled 
mathematically,  and  pointlessly.  A  “longitude  of  perifocus”  is  defined  as  H  +  oj,  which  makes 
no  mathematical  sense  because  D  and  oj  lie  in  different  planes.  (It  is  meaningless  to  add  angles 
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in  different  planes  because  the  sum  retains  no  information  on  the  contribution  from  each  plane, 
and  so  can’t  distinguish  between  different  geometrical  situations.  That  is,  an  orbit  described 
by  Q  =  50°  and  co  =  20°  is  quite  different  to  an  orbit  described  by  fl  =  40°  and  co  =  30°. ) 
Also  a  “mean  longitude”  is  defined  as  D  +  co  +  M0,  a  doubly  meaningless  quantity  because  it 
adds  a  mean  anomaly,  which  is  not  even  an  angle  that  relates  directly  to  the  planet’s  position! 
Despite  their  names,  neither  of  these  two  bizarre  sums  are  angles,  and  they  are  only  ever  used 
as  a  kind  of  secret  code  requiring  decoding  with  a  “key”  f2,  which  is  always  given  and  is  used 
to  extract  co  along  with  M0: 

co  =  long,  of  perifocus  —  D  , 

M0  =  mean  longitude  —  long,  of  perifocus.  (7-1) 

I  won’t  dignify  the  longitude  of  perifocus  and  mean  longitude  with  mathematical  symbols, 
and  I  wish  celestial  navigators  would  stop  using  these  mis-defined  and  mis- named  quantities. 
These  two  parameters  have  no  use  other  than  to  require  decrypting  to  produce  the  actual 
parameters  needed,  co  and  M0  via  (7.1).  Aside  from  their  lack  of  mathematical  meaning,  they 
simplify  nothing,  neither  in  computation  nor  in  understanding  the  subject.  The  longitude  of 
perifocus  and  mean  longitude  are  sometimes  argued  for  on  a  supposed  non-defiuition  of  co  for 
circular  orbits  (which  have  no  perifocus)  or  D  for  equatorial  orbits  (which  have  no  ascending 
node).  But  the  line  of  mathematical  reasoning  that  constructs  an  orbit  does  so  by  starting 
with  a  conic  section  and  referring  it  to  xOP,yOP,  zOP  axes — which  can  certainly  be  done  even 
for  a  circular  orbit  with  its  absence  of  perifocus.  These  axes  then  define  co.  Also,  whilst  an 
equatorial  orbit  has  no  ascending  node,  it  can  still  be  given  a  value  of  D,  in  fact  any  value, 
since  D  serves  only  to  describe  how  the  orbit  is  orientated.  In  such  a  case  the  sum  Q  +  co  is 
well  defined  because  D  and  co  are  in  the  same  plane,  and  increasing  the  choice  of  U  by,  say, 
1°  must  be  offset  by  decreasing  the  choice  of  co  by  1°.  So  only  in  that  case  of  an  equatorial 
orbit  does  Ll  +  co  have  any  meaning — but  it  is  not  needed  there. 

Whether  the  longitude  of  perifocus  and  mean  longitude  have  ever  prevented  a  ship’s  mutiny 
might  never  be  known.  But  those  days  are  long  gone,  and  I  think  the  two  quantities  should 
now  be  relegated  to  the  history  of  orbital  mechanics  [22]. 


8  Relating  Coordinate  Systems 

We  will  need  various  cartesian  coordinate  systems  to  locate  a  planet  in  our  night  sky.  Al¬ 
though  the  following  calculations  can  easily  be  combined  in  a  way  that  enables  some  of  these 
coordinates  to  be  omitted,  it’s  advantageous  to  employ  more  coordinate  systems  than  are 
strictly  necessary,  as  a  way  of  separating  the  calculations  into  manageable  steps. 

First,  orbit-plane  coordinates  centred  on  the  Sun  describe  the  planet’s  motion  most  simply: 
the  planet  moves  right  handed  about  zOP  in  the  x0P  yOP  plane  with  polar  coordinates  r,  8.  The 
perifocus  lies  on  the  positive  xOP  axis.  Take  careful  note  that  orbit-plane  coordinates  are 
planet  dependent,  so  we  write  e.g.  x0PI  for  Jupiter  and  x0PE  for  Earth. 

Next,  Sun-Centred  Inertial  coordinates  SCI  describe  the  unchanging  “global”  Solar  System 
frame,  also  centred  on  the  Sun,  discussed  in  Section  5.  Local  coordinates  such  as  “east-north- 
up”  (ENU)  are  centered  on  the  Earth  observer,  and  are  easily  converted  to  local  bearing 
and  elevation.  We  could  use  just  these  coordinate  sets,  but  will  also  employ  ECI  and  ECEF 
coordinates  as  intermediate  steps. 
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Here  is  a  very  useful  way  to  keep  track  and  make  sense  of  the  various  vectors  and  their 
coordinates  within  these  systems.  First,  remember  that  a  position  is  always  specified  relative 
to  some  point.  To  sight  Jupiter  J  from  Adelaide  A,  we  require  the  position  vector  rJA  of 
Jupiter  relative  to  Adelaide.  The  orbital  parameters  for  Jupiter  will  give  us  the  position  of 
Jupiter  rJS  relative  to  the  Sun  S,  and  the  position  rES  of  Earth’s  centre  E  relative  to  the  Sun. 
Although  it  isn’t  necessary  for  a  calculation  involving  Jupiter,  we  will  carefully  distinguish 
Earth’s  centre  from  the  location  of  Adelaide,  so  will  require  the  position  of  Adelaide  rAE 
relative  to  Earth’s  centre.  This  last  vector  is  really  quite  negligible,  but  we  include  it  to 
show  how  the  various  pieces  of  the  jigsaw  fit  together.  (For  sighting  Jupiter  from  Adelaide, 
the  angular  error  caused  by  omitting  rAE  is  quickly  estimated  as  Earth’s  radius  divided 
by  Jupiter’s  distance  from  the  Sun,  or  approximately  6400  km/(780  million  km),  or  about 
0.0005°.) 

The  above  vectors  are  related  via 


rJA  =  rJS  +  rSE  +  rEA 

=  r.JS  ~  rES  ~  rAE  ■  (8.1) 

Notice  that  no  numbers  are  present  in  (8.1):  it  holds  independently  of  any  coordinate  choice. 
rjA  is  a  proper  vector ,  an  arrow,  but  it  has  coordinates  in  any  coordinate  system  C  that 
we  choose.  These  numbers  form  a  3-elenrent  coordinate  vector  [rjE]c,  which  we  write  as  a 
column  of  numbers,  useful  later  for  multiplying  it  by  a  matrix.  This  notation  was  introduced 
in  Section  2. 

Our  central  task  is  to  relate  the  coordinates  of  rJE  in  any  two  different  cartesian  coordinate 
systems  C  and  C' .  That  is,  given  the  C-coordinates  of  some  vector  v,  what  are  its  C' - 
coordinates  [v  ]e'?  The  following  analysis  is  a  standard  technique  of  linear  algebra,  but  is 
included  here  for  completeness.  The  use  of  just  two  dimensions  here  saves  space,  but  the 
same  results  apply  in  three  dimensions. 

Write  the  unit  basis  vectors  of  C  as  u x ,  uy  and  the  unit  basis  vectors  of  C'  as  uxi ,  uy/ .  The 
following  procedure  demands  each  set  of  basis  vectors  to  be  orthonormal ,  meaning  its  vectors 
have  unit  length  and  are  mutually  orthogonal.  Starting  with 

V  VX'liX  T  VyUy  VXl1XXl  T  VylUyl  ,  (8.2) 


the  relevant  coordinate  vectors  of  v  are 


<2 

_ i 

vux 

IT}] 

1 

i - 

Is 

s 

vv\ 

V  ■  Uy 

5  IWJC 

K'J 

V  ■  Uy 1 

Now  consider 


{VXUX  +  VyUy)  -Uxt 

(yx'U'x  ^ y ' 

Ux  •  UXI  Uy  ■  Uxl 

UX-Uyl  Uy'Uyl 

That  is,  the  matrix  transforms  coordinates  from  C  to  C': 

Me'  =  Me , 


Me'  = 


V  ■  Uyl 


V'Uyl 


(8.3) 


(8.4) 


(8.5) 
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where  (8.4)  shows  that  the  columns  of  are  the  basis  vectors  of  C  expressed  in  C'  coordin¬ 
ates: 


iux}c'  [ uy\c‘ 


(8.6) 


We  will  use  the  idea  of  (8.6)  to  construct  the  necessary  fj,  matrices  ahead  [23].  Equivalently, 
(8.4)  shows  that  the  rows  of  f/f,  are  the  basis  vectors  of  C  expressed  in  C  coordinates.  From 
this,  it’s  not  hard  to  see  that  is  both  the  transpose  and  the  inverse  of  : 


Pc 


(8.7) 


So  a  fi  matrix  is  inverted  by  transposing  it.  These  matrices  allow  a  switch  between  any 
coordinates  that  we  prefer  to  use.  They  can  be  “chained”  together;  for  example,  with  three 
coordinate  systems: 

Ma  =  Pa  M b  =  Pa  P°b  Me  =  Pa  Me .  (8-8) 

showing  that  /i^  =  / 1 The  extension  to  more  coordinate  systems  is  immediate,  so  that 
for  coordinate  systems  A,  B, . . . ,  Y,  Z, 

Pa  =  Pa  Pb  Pc  •  ■  •  Py  ■  (8-9) 


Calculating  a  /x  matrix  presents  a  slight  complication:  the  approach  that  might  be  con¬ 
sidered  as  intuitive  requires  more  calculation  than  is  needed  by  a  particular  less  intuitive 
approach.  Fortunately  we  can  always  convert  each  of  these  approaches  to  the  other,  depend¬ 
ing  on  whether  we  seek  more  intuition  or  more  mathematical  simplicity  (which  equates  to 
computational  speed  here).  As  an  example,  use  (8.6)  as  a  guide  to  calculate  /x°c:,  which 
converts  orbit-plane  coordinates  to  Sun-Centred  Inertial  coordinates: 


Psci  = 


u. 


2/opJSCI 


(8.10) 


Before  we  describe  the  rotations,  it  helps  the  visualisation  to  realise  that  a  proper  vector 
connecting  two  well-defined  points  is  not  tied  to  the  origin  of  any  set  of  coordinates;  you  can 
translate  that  vector  anywhere  at  all  if  that  helps  you  to  visualise  what  we  are  doing.  So  to 
picture  a  rotation  through  some  angle,  imagine  the  vectors  as  moving  over  the  surface  of  a 
sphere  through  that  angle,  all  the  while  being  rigidly  embedded  in  the  local  tangent  plane  to 
the  sphere’s  surface;  this  is  richer  than  simply  imagining  them  to  be  rotating  with  their  tails 
anchored  to  some  point  on  the  rotation  axis.  When  the  rotations  involve  latitude/longitude, 
you  can  also  picture  the  vectors  moving  over  Earth’s  surface  along  small  or  great  circles,  while 
being  rigidly  embedded  in  the  surface’s  tangent  plane;  although  Earth  isn’t  quite  spherical, 
the  commonly  used  definition  of  latitude,  geodetic,  ensures  that  moving  a  vector  in  this  way 
along  a  meridian  between  two  latitudes  will  rotate  it  correctly  through  the  angular  difference 
of  the  latitudes;  for  example,  this  is  why  equation  (8.26)  ahead  is  exact. 

Referring  to  Figure  6,  imagine  taking  copies  of  the  SCI  basis  vectors  ux  ,uy  ,uZsci  and 
rotating  these  copies  through  the  angles  i,Q,u  to  become  the  orbit-plane  set  uXop,uyop,uZop. 
(Remember  that  all  rotations  are  right  handed  unless  otherwise  stated.)  Take  care  to  get 
the  rotation  order  right,  because  two  rotations  cannot  generally  be  swapped.  Picture  the 
sequence  carefully  [24]:  first  rotate  each  copy  of  ux  ,  uysci,  u2gc[  around  zsci  by  U.  then  rotate 
the  results  around  the  rotated  version  of  uXsci  by  i,  and  finally  rotate  the  results  around  the 
doubly  rotated  version  of  uZsa  by  ui,  producing  uXop,uyop,uZop.  Write  this  procedure  as 

■fMzsci’  U2/sci>  'UZsci}  ^  RZs«  *  R(xSCj)  *  R(zscl) 

>  {Uxop’  uyop1  Uzop}  >  (8-11) 
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where  Rzsci  means  “rotate  about  zSCI  (or  uZsrJ  by  0”,  and  (xSCI)  denotes  the  latest-rotated 
copy  of  the  xSCI  basis  vector. 

This  procedure  is  easy  to  visualise,  but  requires  the  carried-along  axes  to  be  calculated 
because  we  are  to  rotate  around  them.  That’s  not  difficult  to  do,  but  much  easier  and  much 
faster  computationally  is  to  realise  that  the  rotations  can  be  done  around  the  original  SCI 
axes/basis  vectors  only:  first  rotate  the  copies  of  uXsci,uy  ,  u%ri  around  zSCI  by  u>,  then  rotate 
the  results  around  xSCI  by  i,  and  finally  rotate  the  results  around  zSCI  by  Q: 

Ir/T  ,  u.j  ,  u ..  I  — R~  — y  RT  — y  R 

L  *^SCI  f/sci  ^SCI  J  ^SCI  ^SCI  "SCI 

~ ^  {u*OP’us/op’Uz0p}‘  (8-12) 


This  last  description  is  perhaps  not  so  intuitively  obvious  as  the  first  one  in  (8.11).  Compare 
(8.11)  and  (8.12):  equation  (8.12)  reverses  the  order  of  the  angles  rotated  through  in  (8.11) 
while  changing  the  sense  of  its  rotation  axes  from  “copies  carried  along”  to  “fixed  originals”. 
It’s  quite  easy  to  show  that  this  reversal  can  always  be  done  for  any  number  of  rotations, 
even  if  the  rotation  axes  are  not  mutually  perpendicular.  This  apparently  nameless  theorem 
simplifies  many  calculations  in  orientational  analysis  but,  surprisingly,  is  not  well  known. 
Reference  [25]  discusses  it  in  detail. 

Rotating  a  proper  vector  about  an  x  axis  though  angle  9  is  accomplished  by  multiplying 
its  coordinate  vector  by  the  Euler  matrix  Ef,  and  similarly  for  rotations  about  the  y  and 
z  axes  [26]: 


1 

0 

0 

11 

cos  0 

0 

sin  6 

,  El  = 

cos  0 

— sin# 

0 

0 

cos  0 

— sin# 

0 

1 

0 

sin  6 

cos  6 

0 

0 

sin  0 

cos  6 

—sin  0 

0 

cos  6 

0 

0 

1 

(8.13) 


Recall  that  we  are  calculating  as  an  example  of  the  method  of  constructing  these  ft  matrices. 
Refer  now  to  (8.10)  and  realise  that  each  column  of  is  one  of  the  OP  basis  vectors  ex¬ 
pressed  in  SCI  coordinates.  So  follow  the  procedure  of  (8.12)  to  construct  the  OP  basis 
vectors,  all  the  while  working  in  SCI  coordinates.  Three  rotations  are  required,  carried  out 
by  three  matrix  multiplications  (note  the  correct  order!): 

o  ■  r°i 

_  ZTi  1  TPl  TT'LJ 

■i  —  -e/3  E1  E3  i  , 

_0_ 

o  ,  (8-14) 

l. 

from  which  matrix  block  multiplication  makes  it  apparent  that 


\u„ 


Jsci 


_  771^2  Tpi  rpt 

—  .C/3  H/i  JL> 


U. 


2/op-Is 


_  rpQ  TT'l  ZUC. 

—  C/Q  C/1  C/Q 


Msci  = 


[U3/0pJsCI  iUz0 


_  ipQ  Tpi  TPU 

—  -C'l  ^3 


10  0 
0  10 
0  0  1 


_  TpQ  Tpi  TPU 

—  ^3  E 1  E3 


(8.15) 


8.1  Combining  the  Various  Coordinate  Systems 

The  bearing  and  elevation  of  Jupiter  in  Adelaide’s  sky  [27]  can  easily  be  derived  from  the 
ENU  coordinates  of  the  position  vector  of  Jupiter  relative  to  Adelaide,  so  we  require  to 
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calculate  [rJA]ENV  by  using  (8.1): 

\rjA  ENU  =  [rjs  ~  rEs  -  rAE\  ENU 


=  Menu  [rJS] op,  ~  iCu  [^S]ope  ~  /Cf  [rAE] 


(8.16) 


where  OPJ  denotes  the  orbit-plane  coordinate  system  of  Jupiter  and  OPE  is  likewise  for  Earth. 
Equation  (4.2)  gives  [j-jsIopj  and  [r-Es]opE>  where  each  of  these  are  the  transpose  of  (x0P,  yOP,  0) 
in  that  equation,  with  xOP,yOP  calculated  in  each  instance  from  the  orbital  parameters  of 
Jupiter  and  Earth  respectively.  The  latitude/longitude  of  Adelaide  suffice  to  give  [t\4e]ecef- 
Here  are  the  details  of  calculating  the  three  coordinate  vectors  in  (8.16).  Jupiter’s  orbital 
elements  are  a,  e,  i,  H,  uj,  M0  at  some  epoch  t0.  Calculate  the  semi-minor  axis  length  b  from  a 
and  e  using  the  expression  in  the  text  just  after  (3.2).  Calculate  Jupiter’s  period  T  from  (3.3). 
Convert  the  epoch  time  t0  to  a  JD,  and  do  likewise  for  the  requested  time  t  to  find  t  — 10. 
Then  evolve  the  mean  anomaly  using  Kepler’s  equation  (4.6),  written  as 

M  =  M0  +  27t  (t  —  tg)/T ,  (8.17) 


remembering  to  replace  the  27t  by  360°  if  you  are  working  in  degrees.  Produce  E  by  solving 
E  —  esinE  =  M  (being  careful  to  use  radians  here).  Now  invoke  (4.2)  to  write 


^OPJ 

a  (cos  E  —  e) 

OsIopj  — 

2/op j 

= 

b  sin  E 

0 

0 

(8.18) 


The  same  calculation  using  Earth’s  orbital  elements  produces  [r’Es]0PE-  A  reasonably  ac¬ 
curate  version  of  the  position  vector  of  a  city  (say,  Adelaide)  relative  to  Earth’s  centre  in 
ECEF  coordinates  is  given  by 


\tae]  ECEF 


cos  A  cos  4> 
cos  A  sin  4>  , 
sin  A 


(8.19) 


where  that  city  has  latitude  A  and  longitude  0,  and  Earth  is  assumed  spherical  with  radius  R. 
This  is  sufficient  for  producing  Jupiter’s  sight  direction  from  Adelaide,  but  a  more  accurate 
expression  is  needed  for  working  with  satellites  that  orbit  Earth.  This  more  accurate  version 
models  Earth  as  an  oblate  spheroid  using  the  WGS-84  set  of  Earth  dimensions.  Write  Earth’s 
spheroidal  radii  as 


equatorial:  a  =  6,378,137  m,  polar:  b  =  6,356,752.3142  m.  (8.20) 


Location  A  lies  at  latitude  A,  longitude  4>,  and  height  h  above  the  WGS-84  spheroid  (that  is, 
h  is  approximate  height  above  mean  sea  level).  Write 


k=\/ 


a 2  cos2  A  +  b2  sin2  A  . 


Then  A  has  a  position  relative  to  Earth’s  centre  of 


[vae] 


ECEF 


(a2 /k  +  h)  cos  A  cos  ^ 
(a2 /k  +  h)  cos  A  sin  0 
(b2/k  +  h)  sin  A 


(8.21) 


(8.22) 
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Figure  7:  Relative  orientations  of  the  ECEF  and  enu  axes.  The  colours  of  the  ecef  axes  here 
match  those  in  Figure  5. 


The  transformation  matrices  in  (8.16),  Menu)  Menu:  an<^  MenuFj  can  each  be  calculated  with 
one  sequence  of  rotations,  but  for  clarity  and  to  avoid  duplication  of  effort,  we  will  break  the 
calculation  of  /r°Nu  and  Fenv  into  several  steps.  These  two  matrices  are  just  the  two  neces¬ 
sary  instances  of  the  generic  matrix  that  converts  orbit-plane  coordinates  for  planet  X 
(meaning  Jupiter  or  Earth)  to  ENU.  Write 

Menu  =  Menu"  M^ef  Meci  mSE*  ,  (8-23) 

and  calculate  each  of  the  transformation  matrices  on  the  right-hand  side  of  (8.23).  Taken  as  a 
group,  these  calculations  convert  the  coordinate  systems  step-by-step  in  the  following  chain: 

Figure  3  Figure  7 

I - 1  I - 1 

OP  - >  SCI  - >  ECI  - >  ECEF  - >  ENU 


Figure  6  Figure  5 


Calculate  /TenuF; 

Menu  =  Iuiecef]enu  )W2/EcefJenu  L^ZiramJENU 
We’ll  find  it  easier  to  calculate  the  transpose  of  MenuF>  which  is  Mecef: 


\u. 


Mecef  = 


MUenJeCEF  ['UJ/enuJeCEF  [UZenuJeCEF 


(8.24) 


(8.25) 


Relative  orientations  of  the  ecef  and  ENU  axes  are  shown  in  Figure  7.  The  ENU  coordinates 
have  their  origin  at  Adelaide,  at  latitude  A,  longitude  cf.  The  sequence  of  rotations  that  takes 
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copies  of  the  ENU  basis  vectors  and  rotates  them  to  becomes  ECEF  basis  vectors  is  (noting 
carefully  the  initial  order) 

\u,,  ,u7  ,ur  }  — >  R~x  —tRf  — >  {u,T  ,  u,n  ,u7  j.  (8.26) 

L  f/ECEF  7  "ECEF  7  -t'ECEF J  If  ECEF  ^ECEF  L  *"ENU  7  i/ENU  7  "ENU J  v  ' 

(Or  imagine  rotating  with  i?zECEF  then  R^y  \  ■  then  apply  the  pseudo-reversing  theorem.) 
Recall  the  earlier  comment  that  even  though  Earth  isn’t  quite  spherical,  the  fact  that  we  are 
using  geodetic  latitude  ensures  that  (8.26)  is  exact. 

Working  in  ECEF  coordinates,  (8.25)  becomes 


,ENTJ  _  J^<f>  R~^ 


Mecef 


^ECEF  2/ECEF 


0  0  1 
1  0  0 
0  1  0 


E*E^ 


0  0  1 
10  0 
0  1  01 


and  it  follows  that  the  sought-after  transformation  matrix  is 


ECEF  _ 

Menu 


i 

0  0  1 
1  0  0 


EX2  Ep 


(8.27) 


(8.28) 


Calculate  Mecef:  Again  work  with  the  transpose.  The  ECEF  axes  are  rotated  from  the 
ECI  axes  by  the  Greenwich  sidereal  angle  7  around  the  z  axis  shared  by  ECEF  and  ECI 
coordinates,  so 


ECEF 

Meci 


u 


ECEF 


ECI 


u. 


Uecef  -I  ECI 


R1 


\  1  0  01 
010 
Lo  0  1 J 


=  El 


and  the  required  transformation  matrix  is 


(8.29) 


MecLe  =  Ep. 


(8.30) 


Calculate  pp :  Again,  calculate  first: 


Msci  = 


u. 


2/eci  Jsci 


(8.31) 


Referring  to  Figure  3,  tilt  Earth  by  rotating  copies  of  the  SCI  basis  vectors  left-handed 
around  xSCI  by  Earth’s  tilt  r  (i.e.  right-handed  by  —  r),  then  rotate  the  resulting  vectors  left- 
handed  around  the  original  zSCI  through  the  precession  angle  p  =  360 °{t  —  t0)/ 25,770  years 
(i.e.  right  handed  by  —p).  The  rotations  are 


\uT  ,uv  ,u7  )— >  RrT  ->{uT  ,uv  .uy  ),  (8.32) 

L  ^SCI  7  y SCI  7  ^SCI  J  ^SCI  ^SCI  L  ^ECI  7  y ECI  7  ^ECI  J  7  v  y 

(or  consider  precession  followed  by  tilt,  and  apply  the  pseudo-reversing  theorem),  so  that 

Msci  =  RZ  Rpcl  [  0  1  0 1  =  Ep  Ep  .  (8.33) 


Now  transpose  to  produce  the  required  matrix: 

PP  =  E{El 


(8.34) 


Calculate  pP'.  This  was  obtained  earlier  in  (8.15),  where  the  in  that  equation  are 

those  for  planet  X. 
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We  can  avoid  some  duplication  of  effort  by  writing  (8.23)  as 


Me™  =  Menu  M sc*  : 

i _ i 

independent 
of  planet 


(8.35) 


and  so  use  (8.28),  (8.30),  and  (8.34)  to  write 

SCI  _  ECEF  ECI  SCI 

Menu  Menu  Mecef  Meci 

=  00  1  E"7  E{  El .  (8.36) 

Also,  from  (8.15)  (with  subscripts  J  and  E  for  Jupiter  and  Earth), 

Msa  =  E^E[J  E“j  ,  /4PiE  =  E%EE[E  E“e  .  (8.37) 


The  sought-after  coordinate  vector  of  Jupiter  relative  to  Adelaide  in  ENU  coordinates  [rj^]END 
can  now  be  calculated  from  (8.16).  Next,  write  the  components  of  [r*,M]ENU  as  ^enu)  Menu;  -^enuj 
and  from  these  extract  Jupiter’s  bearing  /3  (ground  angle  from  north  to  east)  and  elevation  e 
via 

sin  f3  =  xenv/D  ,  cos  (5  =  yENU/D  ,  tan  e  =  zENV/D  ,  (8.38) 

where  D  =  a/^enu2  +  Menu2  • 


9  So  Where  is  Jupiter? 

We  now  combine  everything  with  an  example:  in  which  direction  must  one  look  from  an 
Adelaide  back  yard  to  find  Jupiter  at  9:00  p.m.  local  (daylight  savings)  time  on  22nd  March 
2014? 

Begin  by  collecting  the  orbital  elements  for  Jupiter  and  Earth;  these  can  be  found  on  the 
Internet  [28].  The  elements  used  here  apply  to  the  J2000.0  epoch,  and  I  have  quoted  a  set 
that  uses  the  “longitude  of  perihelion”  and  “mean  longitude”  just  to  show  how  to  deal  with 
these  two  strange  and  unnecessary  parameters.  Also,  aside  from  the  elements  themselves,  to 
save  space  I  will  write  most  numbers  below  to  two  significant  decimal  places,  but  will  use 
more  decimal  places  in  the  calculations. 

Convert  the  requested  date  to  GMT  and  then  to  a  julian  day,  do  the  same  for  the  epoch, 
and  you’ll  find  that  the  time  elapsed  since  the  epoch  is  5193.9375  days:  this  is  t  —  t0  in  (8.17), 
the  duration  for  which  we  require  to  evolve  the  planet’s  position  from  its  initial  position  that 
was  given  in  its  orbital  elements.  Do  this  by  applying  Kepler’s  equation  to  Jupiter  and  Earth 
in  turn,  as  follows. 

The  Calculation  for  Jupiter 

Referring  to  (8.16),  calculate  ^Enu  from  (8.35)  and  [rjg]0Pj  from  Kepler’  theory,  (8.18). 
Jupiter’s  mass  and  orbital  elements  are 

mass  =  1.8986  x  102'  kg, 
a  =  5.203  363  01  AU,  e  =  0.048  392  66, 
i  =  1.305  30°,  n  =  100.556  15°, 
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“long,  of  perihelion”  =  14.753  85°, 
“mean  longitude”  =  34.40438°. 


(9.1) 


(The  semi-major  axis  a  is  given  in  “astronomical  units”,  where  1  AU  is  Earth’s  mean  distance 
from  the  sun,  or  1.495  978  707  x  1011  metres.)  Without  hesitation,  apply  (7.1)  to  decrypt  and 
delete  the  last  two  quantities  in  (9.1),  replacing  them  with  the  meaningful  parameters 


u  =  -85.80°,  M0  =  19.65°.  (9.2) 

We  use  a  value  for  the  gravitational  constant  of  G  =  6.67384  x  10~n  SI  units.  Convert  the 
semi- major  axis  length  a  to  metres,  then  calculate  the  semi-minor  axis  length  b  =  a\/ 1  —  e2  . 
Calculate  Jupiter’s  period  T  from  (3.3),  using  a  solar  mass  of  M  =  1.989  x  1030  kg  and 
Jupiter’s  mass  m  in  (9.1):  the  result  is  T  =  4332.8  days. 

Apply  (8.17)  to  evolve  Jupiter’s  mean  anomaly  at  epoch,  M0,  to  its  value  at  the  requested 
time,  M.  (Don’t  confuse  this  M  with  the  Sun’s  mass.)  The  result  is  M  =  91.20°.  Solve 
E  —  esinE  =  M  for  Jupiter’s  eccentric  anomaly  E  at  the  requested  time:  E  =  93.97°  (use 
radians  in  the  calculation!).  Apply  (8.18)  to  find  Jupiter’s  position  relative  to  the  Sun  in 
Jupiter’s  orbit-plane  coordinates: 


bUslopj 


-0.91 

7.76 

0 


X  1011  m. 


(9.3) 


Now  use  (8.35)-(8.37)  to  calculate  We  need: 


A  =  —34.9°  [Adelaide’s  latitude],  <f>  =  138.60°  [Adelaide’s  longitude], 


7  =  280.46°+  56.71°  [Greenwich  sidereal  angle  (6.2)],  r  =  23.439°  [Earth’s  tilt], 


360°  x  5193.9375  days 
25,770  years  x  365.25  days/year 


0.20°  [Earth’s  precession]. 


These  give 


,,SCi  (8.36) 
Menu 


"  —0.90 

-0.40 

0.17' 

opj  (8-37) 

J  rsci 

0.97 

-0.25 

0.022' 

-0.25 

0.80 

0.55 

0.25 

0.97 

0.0042 

-0.35 

0.45 

-0.82 

-0.023 

0.0017 

1.00 

(9.4) 


(9.5) 


While  not  necessary  to  the  main  calculation,  the  position  of  Jupiter  relative  to  the  Sun  in  SCI 
coordinates  is 

'-2.86' 


[Os]  sci  =  MsciJ  bUs]  op j 


7.27 

0.034 


x  1011  m. 


(9.6) 


This  coordinate  vector  suggests  two  short  checks.  First,  the  relative  smallness  of  the  third 
element  agrees  with  the  fact  that  Jupiter’s  orbit  plane  nearly  coincides  with  that  of  Earth, 
since  Earth’s  orbit  plane  defines  SCI  coordinates.  Second,  the  length  of  bUslsci  is  5.2  AU, 
which  is  approximately  equal  to  the  value  of  a  in  (9.1) — as  expected,  since  Jupiter’s  orbit  has 
a  small  eccentricity. 

What  we  do  need  is  the  matrix  /Ienu  Li  (8.16):  this  is  simply  the  product  of  the  matrices 
in  (9.5)  in  the  order  listed.  We’ve  finished  with  Jupiter,  and  are  halfway  through  the  main 
calculation. 
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The  Calculation  for  Earth 

Repeat  the  above  steps  for  Earth,  using  its  mass  and  orbital  elements 

mass  =  5.9736  x  1024  kg, 
a  =  1.000  00011  AU,  e  =  0.016  710  22, 
i  =  0.000  05°,  n  =  -11.26064°, 

“long,  of  perihelion”  =  102.94719°, 

(9.7) 

“mean  longitude”  =  100.464  35  . 

Note  that  Earth’s  orbital  inclination  i  is  not  exactly  zero  here  owing  to  astronomers  using  a 
mean  plane  in  their  measurements  of  Earth’s  orbit;  this  is  a  higher-order  correction  to  our 
calculation  that  won’t  concern  us.  Astronomers  give  P  a  value  for  the  same  reason.  The 
calculations  in  this  report  would  ordinarily  be  insensitive  to  a  value  of  P  for  Earth,  but  the 
one  tabulated  in  (9.7)  must  be  used  to  decrypt  the  “longitude  of  perihelion”,  using  (7.1). 

By  following  the  same  steps  as  for  Jupiter,  we  arrive  at  the  position  of  Earth  relative  to 
the  Sun  in  Earth’s  orbit-plane  coordinates: 


irEs]  OPE 


0.28 

1.46 

0 


x  1011  m. 


(9.8) 


The  calculation  of  proceeds  in  the  same  way  as  for  Jupiter,  producing 


Menu  = 


-0.18  0.97  0.17 

0.83  0.060  0.55 

0.52  0.24  -0.82 


(9.9) 


The  calculation  for  Earth  is  finished.  Next,  (8.16)  specifies  and  [taeI-ecw  to  be  calcu¬ 

lated. 


Adelaide’s  Position  Relative  to  Earth’s  Centre 

Converting  between  ECEF  coordinates  and  the  east-north-up  coordinates  centred  on  Adelaide 
is  enabled  via  MenuF>  given  by  (8.28): 


ECEF 

Menu 


-0.66  -0.75  0 

-0.43  0.38  0.82 

-0.62  0.54  -0.57 


(9.10) 


The  final  term  necessary  in  (8.16)  is  the  position  of  Adelaide  relative  to  Earth’s  centre,  in 
ECEF  coordinates.  Equation  (8.22)  gives 


[tae\ 


ECEF 


-3.93 

3.46 


X  106m. 


-3.63 


(9.11) 
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Combining  the  Matrices  and  Coordinate  Vectors 

We  now  insert  all  necessary  matrices  and  coordinate  vectors  into  (8.16)  to  arrive  at 


'-1.66' 

x  1011  m  = 

•^ENU 

^JaIenu  — 

6.21 

2/enu 

3.76 

_  ^ENU_ 

Apply  (8.38)  to  xENU,  yENU,  zENU  to  extract  the  direction  of  Jupiter  in  Adelaide’s  sky: 

bearing  /3  =  345.1°,  elevation  e  =  30.34°.  (9.13) 

Jupiter  appears  fairly  low  down,  just  west  of  north.  Our  result  assumes  no  atmosphere,  but 
it  turns  out  that  atmospheric  refraction  increases  this  particular  elevation  by  only  0.03°. 

We  have  assumed  the  orbital  elements  to  be  constant:  U(t)  =  D(t0)  etc.  They  can  be 
evolved  linearly  over  t  —  t0  using  their  small  tabulated  rates  of  increase,  which  I  have  not 
included  in  this  report.  But  this  refinement  turns  out  to  be  negligible  for  the  short  period  of 
time  over  which  we  are  evolving  the  orbit’s  parameters  from  the  J2000.0  epoch. 

The  direction  we  have  calculated  differs  by  0.15°  from  that  returned  by  the  excellent 
Stellarium  software,  which  gives  a  bearing/elevation  of  344.9°/30.31°  for  the  case  of  no  atmo¬ 
sphere.  Stellarium  adds  empirical  refinements  to  its  calculations:  it  applies  standard  formulae 
to  calculate  parameters  that  do  change  slowly  over  long  time  scales,  such  as  the  length  of  the 
sidereal  day.  For  simplicity,  we  have  treated  this  length  as  constant. 

For  extra  precision  and  to  predict  farther  into  the  future,  several  other  factors  can  be  in¬ 
corporated.  The  easiest  is  atmospheric  refraction  [29].  The  travel  time  of  light  from  Jupiter 
to  Earth  can  also  be  considered — although  this  turns  out  to  be  negligible  for  the  calculation 
above.  (That  Earth  turns  appreciably  during  the  light-transit  time  is  not  relevant;  the  im¬ 
portant  point  is  that  Jupiter  doesn’t  move  far  along  its  orbit  in  that  time.)  Then  there  is 
Earth’s  changing  orientation  due  to  its  non-trivial  nutation,  and  the  fact  that  its  day  is  slowly 
lengthening  due  to  the  Moon’s  tidal  drag.  We  could  also  include  the  ecliptic  plane’s  small 
non-inertiality  in  the  definition  of  the  Sun-Centred  Inertial  frame  SCI,  and  add  various  other 
terms  of  ever-decreasing  size  that  relate  to  the  way  the  Solar  System  evolves  over  time. 

One  can  use  the  above  theory  to  build  a  complex  visual  picture  of  our  Solar  System,  by 
predicting  the  positions  of  planets  as  seen  from  any  other  planet:  the  calculation  simply 
replaces  Earth  with  the  new  home  planet.  From  here  it’s  only  a  short  step  to  predicting  the 
direction  of,  say,  Jupiter’s  moon  Io  as  seen  in  the  sky  at  some  specified  point  on  Neptune’s 
moon  Nereid. 


10  The  Sun  and  Moon 

The  Sun’s  position  as  seen  from,  say,  Adelaide  can  easily  be  predicted.  The  relevant  vector 
is  just  the  reverse  of  Adelaide’s  position  with  respect  to  the  Sun:  rSA  =  —rAS,  and  it  can 
be  coordinatised  in  whatever  coordinate  system  is  useful.  In  particular,  calculating  [t\sti]enu 
results  in  times  of  sunrise  and  sunset  when  you  test  trial  times  for  when  the  Sun’s  elevation 
is  zero — although  atmospheric  refraction  is  more  significant  then  and  should  be  included  in 
the  calculation. 
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What  about  predicting  the  Moon’s  position?  Earth,  Moon,  and  Sun  form  a  genuine  3-body 
system  [30],  but  this  only  really  affects  calculations  of  the  Moon’s  position.  Pinpointing  that 
position  accurately  requires  adding  several  empirical  terms  to  the  scheme  described  above. 
These  terms  are  not  described  here  but,  even  so,  we  can  predict  the  position  to  an  accur¬ 
acy  of  a  degree  or  so  by  treating  the  Moon  and  Earth  as  a  2-body  system  using  the  above 
analysis.  (Unfortunately,  because  the  Moon  is  half  a  degree  wide,  the  result  isn’t  accurate 
enough  to  predict  eclipses.)  The  calculation  is  similar  to  that  for  Jupiter,  but  the  3-body 
nature  of  the  Moon’s  orbit  obliges  us  to  use  a  measured  value  of  its  period  rather  than 
use  (3.3)  to  calculate  the  period.  The  period  is  of  course  roughly  one  month,  but  several 
types  of  month  are  defined,  and  to  use  the  above  calculation  of  anomalies  you  must  use  the 
anomalistic  month,  the  time  between  successive  perigees,  of  27.554550  days.  Also,  the  two 
orbital  elements  D  and  uj  are  now  no  longer  even  remotely  constant.  The  Moon’s  ascending 
node  regresses  within  the  ecliptic  plane  (i.e.  moves  left  handed  about  zSCI)  through  360°  every 
18.61  years,  so  that  dhl/dt  =  — 360°/(18.61  years).  If  you  stand  on  this  moving  node  and 
watch  the  Moon’s  orbit  plane,  you’ll  notice  that  the  Moon’s  perigee  rotates  within  that  orbit 
plane  through  360°  every  6.00  years  in  the  same  direction  that  the  Moon  orbits,  resulting  in 
dui/dt  =  360°/ (6.00  years)  [31].  These  two  rates  of  increase  are  easily  used  to  evolve  U  and  ui 
from  the  epoch  to  the  requested  time.  Then  just  apply  (8.15)  as  usual,  using  Q(t)  and  ui(t) 
instead  of  U(t0)  and  uj(t0). 


11  Finding  Stars  and  Nebulae 

Finally,  it  becomes  a  simple  matter  to  use  a  subset  of  the  above  analysis  to  calculate  the 
bearing  and  elevation  of  any  “fixed”  celestial  object  seen  from  any  place  and  time  on  Earth. 
Stellar  positions  are  tabulated  as  the  celestial  longitude  and  latitude  of  their  directions  from 
Earth  with  reference  to  Earth- Centred  Inertial  coordinates:  that  is,  the  sight  direction’s  angle 
east  of  the  First  Point  of  Aries  in  Earth’s  equatorial  plane  xEC1yEC1,  and  its  angle  north  of 
Earth’s  equatorial  plane.  (The  extension  of  Earth’s  equatorial  plane  out  to  a  great  circle  in 
the  sky  produces  the  celestial  equator ,  and  the  vanishing  points  of  Earth’s  rotation  axis  in  the 
sky  are  the  celestial  poles.)  Astronomers  call  an  object’s  celestial  longitude  its  right  ascension, 
and  the  celestial  latitude  its  declination.  We  can  also  calculate  right  ascension  and  declination 
of  the  planets  using  the  analysis  above,  and  so  plot  their  positions  in  a  star  atlas. 


12  How  Far  Ahead  Can  We  Predict? 

If  we  use  the  above  numbers  and  2-body  calculation  to  predict  Jupiter’s  position  for  some 
date  in  AD  3000,  our  result  will  differ  from  that  of  Stellarium  by  a  dozen  degrees.  Clearly 
we  need  more  knowledge  of  the  various  changing  parameters  to  predict  1000  years  into  the 
future.  Evolving  the  orbital  elements  using  their  first  time  derivatives  reduces  this  error  to 
about  7°.  In  fact,  these  orbital  elements  and  their  rates  of  increase  are  determined  partly 
from  observation,  and  partly  from  running  computer  simulations  of  planets’  orbits  and  least- 
squares  fitting  the  elements  and  their  rates  of  increase  to  the  output;  so  from  our  point  of 
view  the  exercise  of  increasing  our  precision  becomes  a  little  artificial. 

At  the  end  of  the  day,  the  hunt  for  more  precision  is  not  a  question  of  requiring  more 
physics.  We  have  all  the  physics  we  need,  but  modelling  an  n-body  system  by  pairs  of  2-body 
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solutions  will  only  ever  be  an  approximation.  We  cannot  hope  to  predict  arbitrarily  far  ahead 
with  such  a  model,  unless  we’re  prepared  to  add  an  arbitrarily  large  number  of  empirical  terms 
that  modify  those  2-body  solutions  to  fit  the  true  n-body  solution  And  this  n-body  solution 
is  only  the  result  of  solving  Newton’s  gravity  numerically  in  some  approximation  to  that  most 
elusive  of  constructs:  the  long-lived  inertial  frame. 


13  Concluding  Remarks 

If  you  have  followed  the  calculations  in  this  report,  you’ll  have  gained  more  than  you  might 
realise.  You  will  have  a  better  understanding  of  Earth’s  place  in  our  Solar  System.  You’ll  be 
able  to  tackle  the  mostly  arcane  books  on  orbital  theory,  which  tend  to  derive  much  of  the 
above  theory — when  they  derive  it  at  all — in  more  convoluted  ways  than  I  have  done  here. 
The  fi  matrices  of  Section  8  form  the  absolute  core  of  orientational  analysis,  so  you  will  have 
a  better  understanding  of  literature  such  as  used  in  modern  3D  movie  making  and  aerospace 
6  degree-of- freedom  modelling  [32].  With  some  computer  graphics  skills  you  could  write  your 
own  planetarium  software.  And  you’ll  gain  an  appreciation  for  the  efforts  of  early  astronomers, 
who  helped  create  our  modern  world  without  the  benefit  of  our  modern  mathematical  and 
computational  tools. 
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r  =  fur  +  riir  .  (14-1) 


UNCLASSIFIED 


27 


DST-Group-TN-1594 


UNCLASSIFIED 


What  is  iir  in  terms  of  ur  and  uel  Write  the  derivative  using  infinitesimal  notation  as 


ur 


ur(t  +  d  t)  —  ur(t ) 
d  t 


(14.2) 


Now  consider  ur  “pinned”  to  the  planet,  and  ask  how  it  evolves  as  the  planet  moves  from  time  t  to  t  +  d t. 
Changing  r  has  no  effect  on  ur,  so  consider  in  the  following  figure  only  what  happens  to  ur  as  #  changes: 


In  a  time  dt,  6  increases  by  d#,  causing  ur  to  turn  through  d#.  Thus  the  increment  ur(t  +  df)  —  ur(t)  has 
length  dt?  and,  in  this  limit  of  infinitesimal  rotation,  is  rotated  90°  from  ur;  so  this  increment  vector  is 
just  d#  ug.  Hence  (14.2)  becomes  ur  =  d#  ug/dt  =  8ug,  so  that  (14.1)  is  r  =  rur  +  r9ug.  Now  repeat  this 
idea:  differentiate  this  last  equation  with  respect  to  time,  and  now  calculate  ug  similarly  to  (14.2)  to  give 
ug  =  —9ur.  We  arrive  at  the  sought-after  expression  for  r  in  terms  of  ur  and  ug,  which  is  equation  (2.7) 
once  again: 

r  =  (r  —  rd2)ur  +  (2r9  +  r6)ug  .  (14.3) 

[8]  For  an  alternative  to  this  vector-calculus  approach  of  producing  (2.8),  try  the  approach  of  writing  Lag¬ 
range’s  equations  in  the  variables  r  and  6,  as  Lagrange’s  approach  doesn’t  rely  on  vectors. 

[9]  F.  van  Diggelen  (2009)  A-GPS:  Assisted  GPS,  GNSS,  and  SBAS,  Artech  House.  See  Chapter  8,  Ephemeris 
Extension,  Long-Term  Orbits. 

[10]  J.  Connor  (2005),  Kepler’s  Witch,  HarperOne.  Page  5  cites  an  essay  by  eminent  science  historian  I.  Bern¬ 
ard  Cohen. 

[11]  You’ll  sometimes  encounter  the  idea  that  a  sort  of  fiducial  idealised  satellite  can  be  envisaged  as  moving 
in  a  circular  orbit  at  constant  speed  about  the  Sun,  with  an  angle  from  perihelion  given  at  any  moment 
by  the  mean  anomaly.  This  picture  is  just  that,  but  has  no  physical  or  mathematical  utility. 

[12]  When  using  trigonometric  functions  to  specify  an  angle  implicitly,  we  must  always  specify  two  pieces  of 
information  about  that  angle:  say,  its  sine  and  cosine,  or  any  one  of  its  sine,  cosine,  or  tangent  along 
with  its  quadrant,  etc.  Most  computer  languages  have  a  function  atan2  that  returns  8  when  given  sin# 
and  cos#  as  arguments  (and  be  sure  always  to  check  the  order  that  these  arguments  must  have).  In 
particular,  it’s  incorrect  to  think  that  (4.11)  can  be  expressed  in  the  form  “#  =  tan-1  ”  for  all  #— 
a  wrong  piece  of  mathematics  that  you’ll  often  find  in  books  and  journal  papers.  The  inverse  tangent 
function  simply  isn’t  defined  that  way;  and  making  such  a  mistake  has  the  effect  of  discarding  one  of  the 
two  required  pieces  of  information.  No  inverse  tangent  function  can  be  defined  in  such  a  way,  because 
any  function  that  returns  #  requires  two  pieces  of  information.  The  correct  expression  is  “#  =  tan-1 

(+  7r  if  #  is  in  quadrant  2  or  3)”. 

[13]  I  won’t  complicate  things  by  saying  that  the  julian  day’s  GMT  differs  from  “modern  GMT”  (Coordinated 
Universal  Time,  UTC)  by  some  leap  seconds,  which  need  to  be  carefully  accounted  for  in  high-precision 
calculations.  Several  other  time  scales  are  also  defined  in  this  subject.  For  more  information,  see  http: 
/ /www .bipm . org/en/bipm- services/time scales. 

[14]  Bizarrely  and  deplorably,  the  terms  julian  day  and  julian  date  have  in  recent  years  sometimes  been 
misappropriated  as  names  for  the  ordinal  date,  which  is  simply  a  number  for  the  current  day  of  the  year, 
from  1  to  365  or  366,  which  resets  to  1  at  the  start  of  each  year. 

[15]  I  emphasise  that  12:00  noon  GMT  on  1st  January  4713  BC  in  the  proleptic  julian  calendar  is  the  same 
moment  as  12:00  noon  GMT  on  24th  November  4714  BC  in  the  gregorian  calendar. 
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Frames  and  coordinates  are  different  things.  A  frame  is  physical:  a  choice  of  viewpoint  visualised  as  a 
lattice  within  which  we  view  and  describe  the  physics  of  a  situation.  Coordinates  are  mathematical:  sets 
of  numbers  that  locate  events  relative  to  the  lattice,  and  used  to  perform  calculations.  A  given  frame  can 
be  described  by  any  number  of  different  coordinate  systems  (e.g.  cartesian  and  polar);  and  conversely,  a 
given  coordinate  system  can  be  used  in  any  number  of  different  frames  (e.g.  cartesian  is  typically  used 
for  SCI  and  ECI). 

The  word  solstice  means,  in  essence,  “Sun  stationary”,  denoting  the  two  moments  in  the  year  when  the 
Sun  reaches  its  most  northerly  and  southerly  points  in  the  sky.  The  word  equinox  means  “equal  night”, 
referring  to  the  fact  that  around  the  equinox  as  the  Sun  crosses  the  celestial  equator,  day  and  night  have 
approximately  equal  durations.  These  two  durations  can  only  be  approximately  equal  because  the  Sun 
doesn’t  necessarily  cross  that  equator  at  a  convenient  civil  time. 

The  First  Point  of  Aries  is  often  called  the  vernal  equinox  or  spring  equinox,  both  terms  that  I  see  as 
misleading  and  should  not  be  used.  A  point  in  the  sky  and  a  moment  in  time  are  two  separate  things, 
and  we  need  terms  for  both:  the  First  Point  of  Aries  is  unambiguously  a  point  in  the  sky,  so  the  equinox 
(“equal  night”)  should  be  reserved  for  a  moment  in  time  (when  the  Sun  crosses  the  celestial  equator). 
“Vernal”  means  green,  referring  to  the  season  of  spring  which  begins  in  the  northern  hemisphere  around 
March,  and  this  certainly  refers  to  time;  however,  what  astronomers  often  call  the  “spring  equinox”  is 
a  moment  in  time  that  occurs  in  the  southern  hemisphere’s  autumn.  No  confusion  will  ever  arise  if  you 
refer  to  the  point  in  the  sky  as  the  First  Point  of  Aries,  and  the  moment  in  time  as  the  March  equinox 
(as  well  as  the  June  solstice,  September  equinox,  and  December  solstice). 

Why  does  Earth  precess?  Earth’s  equatorial  bulge  feels  unbalanced  pulls  in  the  inhomogeneous  gravity 
fields  of  both  Sun  and  Moon.  The  resulting  torque  causes  Earth  to  precess  as  it  spins.  In  contrast  to 
a  spinning  toy  top  that  is  trying  to  be  toppled  by  the  combination  of  gravity  and  the  contact  force  of 
the  table  it  sits  on,  Earth  is  trying  to  be  “righted”  by  Sun  and  Moon,  so  it  precesses  in  the  opposite 
direction  to  that  of  a  spinning  top.  Earth’s  bulge,  along  with  displaced  water  in  its  tides  that  drag  along 
Earth’s  surface  slightly  out  of  phase  with  the  Moon  that  causes  them,  gives  Earth  a  slightly  non-radial 
gravitational  field  that  applies  a  tiny  sideways  force  to  the  Moon,  which  doesn’t  quite  orbit  in  Earth’s 
equatorial  plane.  This  tiny  force  “micro  sling-shots”  the  Moon  away  from  Earth,  resulting  in  the  Moon’s 
distance  from  us  increasing  by  several  centimetres  per  year. 

US  Naval  Observatory  Circular  No.  179  (2005),  freely  available  from  their  web  site. 

Besides  being  used  as  a  measure  of  angle  by  astronomers,  the  sidereal  hour  is  sometimes  said  to  denote  a 
measure  of  time  relative  to  the  stars  instead  of  the  Sun:  24  sidereal  hours  is  said  to  equal  one  sidereal  day, 
or  about  23.9345  civil  clock  hours.  Those  who  insist  on  this  usage  now  have  two  kinds  of  temporal  hour 
and  one  angular  hour.  When  the  angular  hours  are  subdivided  into  minutes  and  seconds,  these  latter  are 
not  even  remotely  related  to  the  commonly  used  minutes  and  seconds  of  arc.  Forty  years  on  from  first 
encountering  these  astronomical  concepts,  I  have  yet  to  find  any  use  for  such  artificial  complexity  in  what 
is  really  a  straightforward  subject. 

What  comes  across  as  a  contrariness  in  the  continued  usage  of  these  two  sums  seems  to  be  part  of  the 
culture  of  modern  astronomy.  I  have  already  mentioned  the  name  “sidereal  time”  for  angle,  along  with 
its  two  definitions  of  “hour”  that  are  often  both  used  in  the  same  sentence.  But  the  list  goes  on.  Some 
(not  all)  planetary  specialists  define  longitude  as  positive  west  of  a  planet’s  prime  meridian,  against  all 
mathematical  and  geophysical  usage.  And  astronomers  who  specialise  in  measuring  radiation  routinely 
plot  a  hot  emitter’s  “power  radiated  per  unit  emitter  surface  area  per  unit  frequency”  versus  not  frequency, 
but  wavelength — which  is  not  mathematically  wrong  per  se,  but  which  does  produce  plots  that  cannot 
be  interpreted  in  the  intuitive  way  of  true  density  plots.  This  practice  requires  the  usual  textbook  maths 
of  blackbody  radiation  to  be  reformulated  (such  as  Wien’s  Law),  to  incorporate  what  is  really  a  misuse 
of  the  concept  of  density. 

Witness  the  field’s  widespread  and  lone  use  of  the  “erg”,  a  unit  of  energy  that  equals  10-7  joules, 
when  it  would  surely  be  far  more  appropriate  to  use  a  large  unit,  such  as  10+24  joules,  when  speaking  of 
exploding  stars.  Similarly,  instead  of  measuring  distances  in  light  years  that  everyone  understands,  some 
astronomers  insist  on  using  the  parsec  (short  for  “parallax-second”),  the  distance  at  which  a  star  has  a 
half-yearly  parallax  of  one  arc  second  when  seen  from  Earth  orbiting  the  Sun.  The  parsec  dates  from  the 
dawn  of  stellar  distance  measurement,  but  its  definition  has  no  modern  use  for  anything  much  farther 
than  a  half  dozen  of  the  nearest  stars;  and  given  that  one  parsec  equals  about  3.26  light  years,  the  unit 
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doesn’t  even  introduce  any  economy  of  numerical  use  as  compared  with  the  light  year.  Quoting  distances 
in  parsecs  nowadays  comes  across  as  an  exercise  in  deliberate  obscurity. 

[23]  The  distinction  between  proper  vectors  and  coordinate  vectors,  along  with  equation  (8.6)  and  the  pseudo- 
reversing  theorem  discussed  after  (8.10),  form  the  heart  of  orientation/rotation  theory.  If  you  digest  what 
I  have  written,  you’ll  be  able  to  perform  and  understand  orientational  calculations  very  easily. 

[24]  Canonical  pictures  of  such  rotation  sequences  sometimes  appear  in  books,  but  I  haven’t  provided  any 
because  you  are  probably  better  off  visualising  my  description  rather  than  attempting  to  work  out  what 
is  going  on  in  those  pictures. 

[25]  Very  few  books  mention  the  pseudo-reversing  theorem,  but  try  proving  it  by  first  twisting  your  hand, 
then  rotating  your  body;  now  return  to  the  start  position  and  rotate  your  body,  then  twist  your  hand. 
Ask  yourself  precisely  which  axes  these  rotations  are  being  done  around,  and  establish  that  any  number 
of  rotations  can  be  “reversed”  in  this  way  when  we  alter  the  sense  of  the  rotation  axes  appropriately. 
You  can  find  this  method  of  reversing  the  rotation  order  explained  in  great  detail  in  D.  Koks  (2012), 
A  Pseudo-Reversing  Theorem  for  Rotation  and  its  Application  to  Orientation  Theory,  DSTO-TR-2675, 
Melbourne,  Vic.,  Defence  Science  and  Technology  Organisation  (Australia). 

[26]  Some  books  change  the  sign  of  6  in  their  Euler  matrices  because  they  define  those  matrices  with  an 
opposite  sense  of  rotation  to  that  used  in  this  report.  Additionally,  notice  that  e.g.  Ef  serves  to  rotate  a 
vector  around  any  generic  x  axis.  For  example,  when  a  vector  v  is  rotated  through  6  around  the  Xa  axis 
of  coordinate  system  A,  the  result  is  v' ,  where 

W]a  =  E{  [v]a  , 

and  when  v  is  rotated  through  6  around  the  xB  axis  of  coordinate  system  B,  the  result  is  v" ,  where 

[v"]B  =  El  [v]B  . 

Hence  we  don’t  write  E l  as  because  there  may  well  be  an  x  axis  present,  and  the  action  of  the 
matrix  is  not  confined  to  rotation  about  only  that  axis. 

[27]  Bearing  and  elevation  are  often  called  “azimuth”  and  “altitude”  respectively  in  astronomy.  Outside  that 
subject,  “azimuth”  doesn’t  necessarily  carry  a  sense  of  being  referred  to  north,  and  of  course  “altitude” 
more  usually  denotes  a  height,  not  an  angle. 

[28]  I  have  used  those  at  http://www.met.rdg.ac.uk/~ross/Astronomy/Planets.html 

[29]  J.  Meeus  (1998),  Astronomical  Algorithms ,  2nd  ed.,  Willmann-Bell.  Equation  (16.4). 

[30]  Earth  and  Moon  are  actually  better  regarded  as  a  double  planet.  Our  Moon  has  exceptional  orbital  and 
physical  features  that  distinguish  it  over  the  other  satellites  in  the  Solar  System,  blurring  its  orbital 
character  to  an  extent  that  means  it  can’t  simply  be  treated  as  Earth’s  satellite.  The  Moon  really  orbits 
the  Sun  in  step  with  Earth,  so  that  calculations  of  its  orbit  must  take  Earth’s  motion  about  the  Sun 
into  consideration.  Strictly  speaking,  a  “keplerian  orbit”  relates  only  to  the  2-body  problem;  but  unlike 
other  satellites  in  the  Solar  System,  the  Moon’s  motion  is  well  and  truly  that  of  a  3-body  problem,  for 
which  no  analytical  solutions  are  known.  Incidentally,  our  Moon  probably  has  more  in  common  with  the 
major  planets  than  the  planet  Pluto  lacks,  which  shows  the  logical  inconsistency  displayed  by  the  mostly 
non-planetary  specialists  of  the  International  Astronomical  Union,  when  a  minority  of  members  decreed 
in  2006  that  we  must  all  no  longer  call  Pluto  a  planet. 

[31]  The  value  of  6.00  years  is  well  known.  Don’t  use  the  figure  of  8.85  years  found  on  many  web  sites.  The 
8.85  years  actually  relates  an  anomalistic  month  to  a  “sidereal  month”,  which  we  are  not  using. 

[32]  For  example,  differentiate  (8.5)  with  respect  to  time  to  produce  what  is  called  the  rotational  derivative 
in  the  reference  that  follows,  a  comparatively  recent  concept  in  aerospace  literature  that  is  actually 
the  aerospace  version  of  the  covariant  derivative  of  vector  and  tensor  calculus.  For  an  introduction  to 
the  rotational  derivative,  see  Zipfel  [3].  That  reference  writes  our  [v\a  =  /r®  [u]s  as  [v]A  =  [T]'4s[u]s. 
(Additionally,  the  use  in  [3]  of  “[ds/dt]'4”  in  an  arbitrary  frame  [not  necessarily  inertial]  must  be  interpreted 
in  that  book’s  notation  as  really  d[s]A/dt  to  be  mathematically  meaningful  [contrast  with  the  comment 
in  [6] ,  which  assumes  inertiality] .  Also,  its  theory  of  the  rotational  derivative  can  be  simplified  notationally, 
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but  is  still  accessible.)  Note  that  the  /r  matrices  in  this  report  are  partly  my  own  notation;  you’ll  find 
all  manner  of  other  notation  in  aerospace  literature.  Most  of  that  literature  doesn’t  distinguish  between 
proper  vectors  (arrows)  and  coordinate  vectors  (arrays  of  numbers) — a  distinction  that  I  see  as  crucial 
for  an  in-depth  understanding  of  calculating  with  vectors  in  multiple  coordinate  systems. 
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